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ABSTRACT 

We investigate the coevolution of galaxies and hosted supermassive black holes throughout the 
history of the Universe by a statistical approach based on the continuity equation and the abundance 
matching technique. Specifically, we present analytical solutions of the continuity equation without 
source term to reconstruct the supermassive black hole (BH) mass function from the AGN luminosity 
functions. Such an approach includes physically-motivated AGN lightcurves tested on independent 
datasets, which describe the evolution of the Eddington ratio and radiative efficiency from slim- to 
thin-disc conditions. We nicely reproduce the local estimates of the BH mass function, the AGN duty 
cycle as a function of mass and redshift, along with the Eddington ratio function and the fraction of 
galaxies with given stellar mass hosting an AGN with given Eddington ratio. We exploit the same 
approach to reconstruct the observed stellar mass function at different redshift from the UV and far-IR 
luminosity functions associated to star formation in galaxies. These results imply that the buildup 
of stars and BHs in galaxies occurs via in-situ processes, with dry mergers playing a marginal role 
at least for stellar masses < 3 x 10 11 M© and BH masses < 10 9 Mq, where the statistical data are 
more secure and less biased by systematic errors. In addition, we develop an improved abundance 
matching technique to link the stellar and BH content of galaxies to the gravitationally dominant 
dark matter component. The resulting relationships constitute a testbed for galaxy evolution models, 
highlighting the complementary role of stellar and AGN feedback in the star formation process. In 
addition, they may be operationally implemented in numerical simulations to populate dark matter 
halos or to gauge subgrid physics. Moreover, they may be exploited to investigate the galaxy/AGN 
clustering as a function of redshift, mass and/or luminosity. In fact, the clustering properties of BHs 
and galaxies are found to be in full agreement with current observations, so further validating our 
results from the continuity equation. Finally, our analysis highlights that (i) the fraction of AGNs 
observed in slim-disc regime, where anyway most of the BH mass is accreted, increases with redshift; 
(ii) already at z > 6 substantial dust amount must have formed over timescales < 10 8 yr in strongly 
starforming galaxies, making these sources well within the reach of ALMA surveys in (sub-)millimeter 
bands. 

Subject headings: black hole physics - galaxies: formation - galaxies: evolution - methods: analytical 
- quasars: supermassive black holes 


1. INTRODUCTION 

Kinematic and photometric observations of the very central regions in local, massive early-type galaxies strongly 
support the almost ubiquitous presence of black holes (BHs) with masses Mbh > 10 6 M© (Dressier 1989; Kormendy 
& Richstone 1995; Magorrian et al. 1998; for a recent review Kormendy & Ho 2013). Their formation and evolution 
is a major problem in astrophysics and physical cosmology. 

The correlations between the central BH mass and galaxy properties such as the mass in old stars (Kormendy & 
Richstone 1995; Magorrian et al. 1998; Marconi & Hunt 2003; Haring & Rix 2004; McLure & Dunlop 2004; Ferrarese 
& Ford 2005; Graham 2007; Sani et al. 2011; Beifiori et al. 2012; McConnell & Ma 2013; Kormendy & Ho 2013), 
the velocity dispersion (Ferrarese & Merritt 2000; Gebliardt et al. 2000; Tremaine et al. 2002; Giiltekin et al. 2009; 
McConnell & Ma 2013; Kormendy & Ho 2013; Ho & Kim 2014), and the inner light distribution (Graham et al. 2001; 
Lauer et al. 2007; Graham & Driver 2007; Kormendy & Bender 2009) impose strong ties between the formation and 
evolution of the BH and that of the old stellar population in the host galaxy (Silk & Rees 1998; Fabian 1999; King 
2005; for a recent review see King 2014). 

A central role in this evolution is played by the way dark matter (DM) halos and associated baryons assemble. So 
far it has been quite popular, e.g. in most semi-analytic models, to elicit merging as the leading process; as to the 
baryons, ‘wet’ and ‘dry’ mergers or a mixture of the two kinds have been often implemented (for a recent review see 
Somerville & Dave 2015). On the other hand, detailed analyses of DM halo assembly indicate a two-stage process: an 
early fast collapse during which the central regions reach rapidly a dynamical quasi-equilibrium, followed by a slow 
accretion that mainly affects the halo outskirts (e.g., Zhao et al. 2003; Wang et al. 2011; Lapi & Cavaliere 2011). Thus 
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one is led to consider the rapid starformation episodes in the central regions during the fast collapse as the leading 
processes in galaxy formation (e.g., Lapi et al. 2011, 2014; Cai et al. 2013). Plainly, the main difference between 
merging and fast collapse models relates to the amount of stars formed in-situ (e.g., Moster et al. 2013). 

While N —body simulations of DM halo formation and evolution are nowadays quite robust (though details of their 
results are not yet fully understood), the outcomes of lrydrodynamical simulations including star formation and central 
BH accretion are found to feature large variance (Scannapieco et al. 2012; Frenk & White 2012). This is expected 
since most of the relevant processes involving baryons such as cooling, gravitational instabilities, angular momentum 
dissipation, star formation and supermassive BH accretion occur on spatial and temporal scales well below the current 
resolution. 

On the other hand, observations of active galactic nuclei (AGNs) and galaxies at different stages of their evolution 
have spectacularly increased in the last decade at many wavelengths. In particular, the AGN luminosity function is 
rather well assessed up to 2 ~ 6 though with different uncertainties in the X-ray (Ueda et al. 2014; Fiore et al. 2012; 
Buchner et al. 2015; Aird et al. 2010, 2015), UV/optical (Richards et al. 2006; Croorn et al. 2009; Masters et al. 
2012; Ross et al. 2012; Fan et al. 2006; Jiang et al. 2009; Willott et al. 2010a), and IR bands (Richards et al. 2006; 
Fu et al. 2010; Assef et al. 2011; Ross et al. 2012); these allow to infer the BH accretion rate functions at various 
redshifts. In addition, luminosity functions of galaxies are now available up to z ~ 10 in the ultraviolet (UV; Bouwens 
et al. 2015; Finkelstein et al. 2014, Weisz et al. 2014; Cucciati et al. 2012; Oesch et al. 2010; Reddy & Steidel 2009; 
Wyder et al. 2005) and up to 2 ~ 4 in the far-infrared (FIR) band (Lapi et al. 2011; Gruppioni et al. 2013; Magnelli 
et al. 2013); these allow to infer the star formation rate (SFR) function at various redshifts. 

As for galaxies selected by their mid- and far-IR emission, the distribution function of the luminosity associated 
to the formation of massive stars shows that at z < 4 the number density of galaxies endowed with star formation 
rates M* > 10 2 yr -1 is N (log M„) > 10 -3 Mpc -3 . The density is still significant, V(logM*) > 10 -5 Mpc -3 , for 
M„ w 10 3 M q yr -1 . On the other hand, the UV selection elicits galaxies forming stars at much lower rates M* < 30 M 0 
yr -1 up to 2 < 10. The complementarity between the two selections is ascribed to the increasing amount of dust in 
galaxies with larger SFRs (Steidel et al. 2001; Mao et al. 2007; Bouwens et al. 2013, 2015; Fan et al. 2014; Cai et 
al. 2014; Heinis et al. 2014). From deep, high resolution surveys with ALMA at (sub-)mm wavelengths there have 
been hints of possible source blending at fluxes SWo/im > 10 mJy (Karim et al. 2013; Ono et al. 2014; Simpson et al. 
2015b). On the other hand, observations at high spatial resolution of sub-nun selected, high redshift galaxies with the 
SMA and follow-ups at radio wavelengths with the VLA show that z < 6 galaxies exhibiting M* ss a few 10 3 Mq yr -1 
have a number density N ~ 10 -6 Mpc -3 (Barger et al. 2012, 2014), fully in agreement with the results of Lapi et al. 
(2011) and Gruppioni et al. (2013) based on Herschel (single dish) surveys. 

Studies on individual galaxies show that several sub-nnn galaxies at high redshift exhibit M* > 10 3 Mq yr -1 
concentrated on scales < 10 kpc (e.g., Finkelstein et al. 2014; Neri et al. 2014; Rawle et al. 2014; Riechers et al. 
2014; Ikarashi et al. 2014; Simpson et al. 2015a; Scoville et al. 2014). Size ranging from a few to several kpc of 
typical high— z strongly star forming galaxies has been confirmed by observations of many gravitational lensed objects 
(e.g., Negrello et al. 2014). In addition, high spatial resolution observations around optically selected quasars put in 
evidence that a non negligible fraction of host galaxies exhibits M* > 10 3 M 0 yr -1 (Omont et al. 2001, 2003; Carilli 
et al. 2001; Pricldey et al. 2003; Wang et al. 2008; Bonfield et al. 2011; Mor et al. 2012). 

The clustering properties of luminous sub-nun selected galaxies (Webb et al. 2003; Blain et al. 2004; Weiss et al. 
2009; Hickox et al. 2012; Bianchini et al. 2015) indicate that they are hosted by large halos with masses Mh > several 
10 12 Mq and that the star formation timescale is around ~ 0.5 — 1 Gyr. 

The statistics on the presence of AGNs along the various stages of galaxy assembling casts light on the possible 
reciprocal influence between star formation and BH accretion (for a recent review, see Heckman & Best 2014 and 
references therein), although the fine interpretation of the data is still debated. On one side, some authors suggest 
that star formation and BH accretion are strongly coupled via feedback processes, while others support the view that 
the two processes are only loosely related and that the final relationships among BH mass and galaxy properties are 
built up along the entire Hubble time with a relevant role of dry merging processes. 

Most recently, Lapi et al. (2014) have shown that the wealth of data at z > 1 strongly support the view that galaxies 
with final stellar mass Af+ > 10 11 Mq proceed with their star formation at an almost constant rate over ~ 0.5 — 1 
Gyr, within a dusty interstellar medium (ISM). At the same time several physical mechanisms related to the star 
formation, such as gravitational instabilities in bars or dynamical friction among clouds of starforming gas or radiation 
drag (Norman & Scoville 1988; Shlosman et al. 1989, 1990; Shlosman & Noguchi 1993; Hernquist & Mihos 1995; 
Noguchi 1999; Umemura 2001; Kawakatu & Umemura 2002; Kawakatu et al. 2003; Thompson et al. 2005; Bournaud 
et al. 2007, 2011; Hopkins & Quataert 2010, 2011), can make a fraction of the ISM loose angular momentum and flow 
into a reservoir around the seed BH. The accretion from the reservoir to the BH can be as large as 30 — 50 times the 
Eddington rate, leading to slim-disc conditions (Abramowicz et al. 1988; Watarai et al. 2001; Blandford & Begelman 
2004; Li 2012; Begelman 2012; Madau et al. 2014; Volonteri & Silk 2014), with an Eddington ratio A < 4 and an 
average radiative efficiency e < 0.1. This results in an exponential increase of the BH mass and of the AGN luminosity, 
with an e—folding timescale r e f ranging from a few to several 10' years. Eventually, the AGNs at its maximum power 
can effectively transfer energy and momentum to the ISM, removing a large portion of it from the central regions and 
so quenching the star formation in the host. The reservoir around the BH is no more fed by additional gas, so that 
even the accretion and the nuclear activity come to an end. 
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More in general, we can implement lightcurves for the luminosity associated to the star formation and to the BH 
accretion in a continuity equation approach. In the context of quasar statistics, the continuity equation has been 
introduced by Cavaliere et al. (1971) to explore the optical quasar luminosity evolution and its possible relation with 
the radiosource evolution. Soltan (1982) and Choksi & Turner (1992) exploited the mass-energy conservation to derive 
an estimate of the present mass density in inactive BHs. The extension and the derivation of the BH mass function 
has been pioneered by Small & Blandford (1992), who first attempted to connect the present day BH mass function 
to the AGNs luminosity evolution. A simplified version in terms of mass-energy conservation has been used by Salucci 
et al. (1999), who have shown that the distribution of the mass accreted onto central BHs during the AGNs activity 
well matches the mass function of local inactive BHs. A detailed discussion of the continuity equation in the quasar 
and central supermassive BH context has been presented by Yu & Lu (2004, 2008). In the last decade the continuity 
equation has been widely used, though the lightcurve of the AGN, one of the fundamental ingredients, was largely 
based on assumptions (e.g., Marconi et al. 2004; Shankar et al. 2009, 2013; Merloni & Heinz 2008; Wang et al. 2009; 
Cao 2010). Results on the BH mass function through the continuity equation have been reviewed by Kelly & Merloni 
(2012) and Shankar (2013). 

We will also implement the continuity equation for the stellar content of galaxies. This has become possible because 
the UV surveys for Lyman Break Galaxies, the wide surveys HerMES and H-ATLAS obtained with the Herschel space 
observatory made possible to reconstruct the star formation rate function in the Universe out to 2 < 6 for SFRs 
M* ~ 10 — 1000 Mq yr^ 1 . Therefore, we can exploit the continuity equation, in an analogous manner as routinely 
done for the AGNs; the BH mass is replaced by the mass in stars and the bolometric luminosity due to accretion is 
replaced by the luminosity generated by the formation of young, massive stars. 

As for the stellar mass function, it is inferred by exploiting the observed luminosity function in the wavelength range 
of the SED dominated by the emission from older, less massive stars. The passage from stellar luminosity to mass 
is plagued by several problems, which result in uncertainties of order of a factor of 2, increasing for young, dusty 
galaxies (e.g., Cappellari et al. 2013; Conroy et al. 2014). Therefore the mass estimate is more robust for galaxies 
with quite low star formation and/or in passive evolution. All in all, the stellar mass function of galaxies is much 
easier to estimate, and hence much better known, than the BH mass function, particularly at high redshift. Reliable 
stellar mass functions are available both for local and redshift up to z ~ 6 galaxy samples (e.g., Bernardi et al. 2013; 
Maraston et al. 2013; Ilbert et al. 2013; Santini et al. 2012; Stark et al. 2009; Gonzalez et al. 2011; Duncan et al. 
2014). The comparison between the observed stellar mass function and the results from the continuity equation sheds 
light on the relative contribution of dry merging and of in-situ star formation. In the present paper we will solve the 
continuity equation for AGN and for the stellar component after inserting the corresponding lightcurves derived from 
the data analysis of Lapi et al. (2011, 2014, see above). 

Once the stellar and the BH mass functions at different redshifts are known, these can also be compared with the 
abundance of DM halos to obtain interesting relationships between halo mass and galaxy/BH properties. Such a 
technique, dubbed abundance matching, has been exploited by several authors (e.g., Vale & Ostriker 2004; Shankar 
et al. 2006; Moster et al. 2010, 2013; Beliroozi et al. 2013; Behroozi & Silk 2015). In this paper, the technique is 
refined and used in connection with the outcomes of the continuity equation to tackle the following open issues in 
galaxy formation and evolution: 

• Is the BH mass function reflecting the past AGN activity? What was the role of merging in shaping it? (see 
Sect. 2.1 and Appendix B) 

• How does the BH duty cycle evolve? What can we infer on the radiative efficiency and on the Eddington ratio 
of active BHs? (see Sect. 2.1.4) 

• Is there any correlation between the central BH mass and the halo mass and how does it evolve with time? (see 
Sect. 3.1.1) 

• Which is the relationship between the AGN bolometric luminosity and the host halo mass? Can we use this 
relationship with the duty cycle to produce large simulated AGN catalogs? (see Sect. 3.1.1) 

• Which are the bias properties of AGNs? Do they strongly depend on luminosity and redshift? (see Sect. 3.1.1) 

• Can the evolution of the stellar mass function be derived through the continuity equation as in the case of the 
BH mass function, by replacing the accretion rate with the SFR? Does dry merging play a major role in shaping 
the stellar mass function of galaxies? What is the role of the dust in the star formation process within galaxies? 
(see Sect. 2.2, App. B and C) 

• Which is the relationship between the SFR, the stellar mass of the galaxies, and the mass of the host halo? 
Does the star formation efficiency (i.e., the fraction of baryons going into stars) evolve with cosmic time? (see 
Sect. 3.1.2) 

• Which is the relationship between the bolometric luminosity of galaxies due to star formation and the host halo 
mass? Can we use this relationship with the stellar duty cycle to produce large simulated catalogs of star forming 
galaxies? (see Sect. 3.1.2) 
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• Which are the bias properties of star forming and passively evolving galaxies? (see Sect. 3.1.2) 

• How does the specific star formation rate evolve with redsliift and stellar mass? (see Sect. 3.1.3) 

• Which is the relationship between the BH mass and the stellar mass at the end of the star formation and BH 
mass accretion epoch? Does it evolve with time? (see Sect. 3.1.4) 

• How and to what extent can we extrapolate the relationships for both galaxies and hosted AGNs to higher, yet 
unexplored, redshift? (see Sects. 3 and 4) 

To answer these questions we have organized the paper as follows. In Sect.[2]we present the statistical data concerning 
AGN and starforming galaxies, introduce and motivate the corresponding lightcurves, and solve the continuity equation 
to derive the BH and stellar mass functions at different redshifts. In Sect. [3] we exploit the abundance matching 
technique to infer relationships among the properties of the BH, stellar, and dark matter component in galaxies. In 
Sect. |4] we discuss and summarize our findings. 

Throughout this work we adopt the standard flat concordance cosmology (Planck Collaboration 2014, 2015) with 
round parameter values: matter density Dm = 0.3, baryon density fib = 0.05, Hubble constant Hq = 100 h km s -1 
Mpc -1 with h = 0.7, and mass variance er 8 = 0.8 on a scale of 8h _1 Mpc. Stellar masses and luminosities (or SFRs) 
of galaxies are evaluated assuming the Chabrier’s (2003) IMF. 

2. CONTINUITY EQUATION 

Given an evolving population of astrophysical sources, we aim at linking the luminosity function N(L, t) tracing 
a generic form of baryonic accretion (like that leading to the growth of the central BH or the stellar content in the 
host galaxies) to the corresponding final mass function N(M,t). To this purpose, we exploit the standard continuity 
equation approach (e.g., Small & Blandford 1992; Yu & Lu 2004), in the integral formulation 

C°° dr 

N(L,t) = J ' AM [d t N(M,t ) - S(M,t)] £ -^±(L\M,t) ; (1) 

here r is the time elapsed since the triggering of the activity (the internal clock, different from the cosmological time 
t), and dr /AL is the time spent by the object with final mass M in the luminosity range [L, L + AL\ given a lightcurve 

L(t\M, t); the sum allows for multiple solution t* of the equation L = L(t\M, t ). In addition, S(M, t ) is a source term 

due to ‘dry’ merging (i.e., adding the whole mass content in stars or BHs of merging objects without contributing 
significantly to star formation or BH accretion). In solving Eq. (JTJ we shall set the latter to zero, and investigate the 
impact of dry merging in the dedicated App. B. Note that by integrating Eq. dT]) in df from 0 to the present time to, 
one recovers Eq. (18) of Yu & Lu (2004). 

If the timescales t* (that encase the mass-to-energy conversion efficiency) are constant in redshift and luminosity, 
then a generalized Soltan (1982) argument concerning the equivalence between the integrated luminosity density and 
the local, final mass density can be straightforwardly recovered from Eq. © without source term by multiplying by L 
and integrating it over L and t 

f'to pOO poo poo l poo 

/ At / ALLN(L,t) = / AMN(M,t) / ALL V —i = const x / AMMN(M,t) , (2) 

Jo Jo Jo Jo i dT Jo 

where the last equivalence holds since X]?:/ dr,; L = const xM. Specifically, for the BH population the constant is 
equal to e c 2 /(1 — e) in terms of an average radiative efficiency e ~ 0.1. We shall see that an analogous expression holds 
for the stellar component in galaxies. 

More in general, Eq. © constitutes an integro-differential equation in the unknown function N(M,t), that can be 
solved once the input luminosity function N(L,t) and the lightcurve L{t\M, t) have been specified. Specifically, we 
shall use it to derive the mass function of the supermassive BH and stellar component in galaxies throughout the 
history of the Universe. Remarkably, we shall see that Eq. © can be solved in closed analytic form under quite 
general assumptions on the lightcurve. 

2.0.1. Connection with standard approaches for BHs 

It is useful to show the connection of Eq. © with the standard, differential form of the continuity equation for 
the evolution of the BH mass function as pioneered by Small & Blandford (1992) and then used in diverse contexts 
by many authors (e.g., Salucci et al. 1999; Yu & Tremaine 2002; Marconi et al. 2004; Shankar et al. 2004, 2009; 
Merloni & Heinz 2008; Wang et al. 2009; Cao 2010). Following Small & Blandford (1992), BHs are assumed to 
grow in a single accretion episode, emitting at a constant fraction A = Taon/TeiM of their Eddington luminosity 
Le<m = Mbrc 2 /tEdd ~ 1-38 x 10 38 Mbh/A/q erg s _1 in terms of the Eddington time tEdd ~ 4 x 10 s yr. The resulting 
lightcurve can be written as 

_ ^ ^BH £ e (r-Tlif e )/Tef 
^Edd 


Tagn(t|M B h, t) 


T < 7iif e ; 


(3) 
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here Mbh is the final BH mass, ruf e = f dr is the total duration of the luminous accretion phase, and r e f = e tEdd/A (1 — 
e) is the e—folding time in terms of the mass-energy conversion efficiency e. Then one has 

-jr - = T ef ©if [AaGN < LaGn(M B h)] , (4) 

d-t^AGN dAGN 

where the Heaviside step function O h(-) specifies that a BH with final mass Mbh cannot have shone at luminosity ex¬ 
ceeding Tagn(Mbh) = XMbhc 2 /t-Edd- Equivalently, only BHs with final masses exceeding Mbh(Lagn) = Lt^dd/Xc 2 
can have attained a luminosity Lagn and so can contribute to the integral on the right hand side of Eq. ©■ Hence 


such an equation can be written as 

nOO 

Lagn IV(Aagn, t) = / dM B H r e f . (5) 

J Mbh(^agn) 

Differentiating both sides with respect to L and rearranging terms yields 

dtN(M B H,t) - 9m bh [-^AGN A r (TAGN,t)]|L AGN (M BH ) = ^(A^bh,^) ■ (6) 

Tef 

Now one can formally write that 

dilfRH 

IV (Lagn, i) = IV(Mbh, t) 77 -(5agn(Mbh,Q (Lagn} = <5agn(Mbh, t) Lagn (7) 

d-OAGN 

in terms of the BH duty cycle Jagn(Mbh,I) = T\a e {M B Yi,t)/t < 1. Since by definition (Lagn) = c(Mbh)c 2 /(1 — e), 
one finally obtains the continuity equation in the form 

d t N(MBu,t) + 9m bh [(Mbh) IV(M B h,£)] = ; ( 8 ) 


the underlying rationale is that, although individual BHs turn on and off, the evolution of the BH population depends 
only on the mean accretion rate (Mbh}- 


2.1. The BH mass function 

We now solve Eq. m to compute the BH mass function at different redshifts. 


2.1.1. BHs: luminosity function 

Our basic input is constituted by the bolometric AGN luminosity functions, that we build up as follows. We start 
from the AGN luminosity functions at different redshifts observed in the optical band by Richards et al. (2006), Groom 
et al. (2009), Masters et al. (2012), Ross et al. (2012), Fan et al. (2006), Jiang et al. (2009), Willott et al. (2010a), 
and in the hard X-ray band by Ueda et al. (2014), Fiore et al. (2012), Buchner et al. (2015), Aird et al. (2015). 

Then we convert the optical and X-ray luminosities to bolometric ones by using the Hopkins et al. (2007) corrections^. 
Note that in the literature several optical and X-ray bolometric corrections have been proposed (see Marconi et al. 
2004; Hopkins et al. 2007; Shen et al. 2011; Lusso et al. 2012; Runnoe et al. 2012). Those by Marconi et al. (2004) 
and Lusso et al. (2012) are somewhat smaller by < 40% in the optical and by < 30% in the hard X-ray band with 
respect to Hopkins et al. (2007). In fact, since bolometric corrections are intrinsically uncertain by a factor ~ 2 
(e.g., Vasudevan & Fabian 2007; Lusso et al. 2012; Hao et al. 2014), these systematic differences between various 
determinations are not relevant. We shall show in Sect. l2.lHl that our results on the BH mass function are marginally 
affected by bolometric corrections. In addition, we correct the number density for the fraction of obscured (including 
Compton thick) objects as prescribed by Hopkins et al. (2007) for the optical data and according to Ueda et al. (2014; 
see also Ueda et al. 2003) for the hard X-ray data. We stress that both the bolometric and the obscuration correction 
are rather uncertain, with the former affecting the luminosity function mostly at the bright end, and the latter mostly 
at the faint end. 

Given the non-homogeneous nature and the diverse systematics affecting the datasets exploited to build up the 
bolometric luminosity functions, a formal minimum y 2 —fit is not warranted. We have instead worked out an analytic 
expression providing a sensible rendition of the data in the relevant range of luminosity and redshift. For this purpose, 
we use a modified Schechter function with evolving characteristic luminosity and slopes. The luminosity function in 
logarithmic bins N(log Lagn) = IV(Lagn) -^agn ln(10) writes 


IV (log Lagn, 2 ) 




Lagn 

L c (z) 


l-a(z) 


exp 


Lagn 

Lc(z)\ J 


(9) 


7 Most of the optical data are given in terms of magnitude M 1450 at 1450 A. First, we convert them to B— band (4400 A) using the 
relation Mg = M 1450 — 0.48 (Fan et al. 2001), then we pass to B— band luminosities in solar units logZ/g/Z/g,© = —0.4 (Mg — 5.48), and 
finally we go to bolometric luminosities in solar units after Lagn/^o = k g Lb/Lb,q x Lb,q/Lq. For this last step we recall that the 
B— band luminosity of the Sun Lg q « 2.13 x 10 33 erg s —1 « L©/2 is about half its bolometric one L© « 3.9 x 10 33 erg s —1 . In some 
other instances the original data are expressed in terms of a z = 2 K- corrected i— band magnitude Mi(z = 2). We adopt the relation with 
the 1450 magnitude M 1450 = Mf(z = 2) + 1.486 (Richards et al. 2006) and then convert to bolometric as above. 
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The normalization log$(z), the characteristic luminosity logL c (.z), and the characteristic slopes a(z) and u)(z) evolve 
with redshift according to the same parametrization 

p(z) =Po + k p iX + k P 2 X 2 + k P 3 X 3 (10) 

with x = log[(l + z)/(l + zo)] and 20 = 0.1. The parameter values are reported in Table 1. The functional form adopted 
here is similar to the widely-used double powerlaw shape (e.g., Ueda et al. 2014; Aird et al. 2015), but with a smoother 
transition between the faint and bright end slopes; all in all, it provides a data representation of comparable quality. 
In fact, we stress that the results of the continuity equation approach are insensitive to the specific parameterization 
adopted for the luminosity function and its evolution, provided that the quality in the rendition of the data be similar 
to ours. For example, in Sect. 12.1.41 we shall show explicitly that our results on the BH mass function are marginally 
affected when using a double powerlaw shape in place of a modified Scliechter to represent the AGN luminosity 
functions. 

In Fig. |T] we illustrate the bolometric AGN luminosity function at various redshifts, including both our data col¬ 
lection and our analytic parameterization of Eq. ©, with an estimate of the associated ler uncertainty; the z = 10 
extrapolation is also shown for illustration. In the inset we plot the evolution with redshift of the AGN luminosity 
density, computed as 

Pl AG¥S {z) = j dlogLAGN A r (logLAGN, z) Tagn (11) 

and the contribution to the total by specific luminosity ranges. 

2.1.2. BHs: lightcurve 

As a further input to the continuity equation, we adopt the following lightcurve (Yu & Lu 2004) 

^agnM-Mbh^) = ^0 -Mbh.p c 2 /tEdd e^ r_Tp ^ Tef 0 < r < r P 

= AoM B H,pc 2 /^Edde - ^ T-Tp ^ TD t p < t < t p + £td (12) 

= 0 T > Tp + £t d 

This includes two phases: an early one up to a peak time Tp when the BH grows exponentially with a timescale r e f to 
a mass M bh,p and emits with an Eddington ratio Ao until it reaches a peak luminosity Ao -Mbh.p c 2 /tpdd', a late phase 
when the luminosity declines exponentially on a timescale Tp, up to a time 7p + £td when it shuts off. With Ao, eo, 
we denote the average Eddington ratio and radiative efficiency during the early, ascending phase. The e—folding time 
associated to them is r e f = eo ^Edd/Ao (1 — eo). 

The lightcurve in Eq. m has been set in Lapi et al. (2014) in order to comply with the constraints imposed by 
large observational datasets concerning: 

• the fraction of X-ray detected AGNs in FIR/K-band selected host galaxies (e.g., Alexander et al. 2005; Mullaney 
et al. 2012a; Wang et al. 2013a; Johnson et al. 2013); 

• the fraction of FIR detected galaxies in X-ray AGNs (e.g., Page et al. 2012; Mullaney et al. 2012b; Rosario et 
al. 2012) and optically selected quasars (e.g., Mor et al. 2012; Wang et al. 2013b; Willott et al. 2015); 

• related statistics via stacking of undetected sources (e.g., Basu-Zych et al. 2013). 

The same authors also physically interpreted the lightcurve according to a specific BH-galaxy coevolution scenario. 

In a nutshell, the scenario envisages that the early growth of the BH occurs in an interstellar medium rich in gas and 
strongly dust-enshrouded (Lapi et al. 2014; also Chen et al. 2015). The BH accretes in a demand-limited fashion with 
values of Eddington ratios A appreciably greater than unity, though the radiative efficiency e may keep to low values 
because slim-disc conditions develop. Since the BH mass is still small, the nuclear luminosity, though appreciably 
super-Eddington, is much lower than that of the starforming host galaxy, and the whole system behaves as a sub-mm 
bright galaxy with an X-ray nucleus. On the other hand, close to the peak of the AGN lightcurve, the BH mass 
has grown to large values, and the nuclear emission becomes comparable or even overwhelms that of the surrounding 
galaxy. Strong winds from the nucleus remove gas and dust from the ambient medium stopping the star formation 
in the host, while the whole system shines as an optical quasar. If residual gas mass is still present in the central 
regions, it can be accreted in a supply-driven fashion so originating the declining part of the lightcurve; this phase 
corresponds to the onset of the standard thin disk accretion, which yields the observed SEDs of UV/optically-selected 
type-1 AGNs (Elvis et al. 1994; Hao et al. 2014). Actually, the data concerning the fraction of starforming galaxies 
in optically-selected quasar samples suggest such a descending phase to be present only for luminous objects, while 
in low-luminosity ones tiny residual mass is present and the AGN fades more drastically after the peak. When the 
accreting gas mass ends, the BH becomes silent, while the stellar populations in the galaxy evolve passively. For the 
most massive objects, the outcome will be a local elliptical-type galaxy with a central supermassive BH relic. 

All in all, we set the fiducial values of the parameters describing the BH lightcurve on the basis of the Lapi et al. 
(2014) analysis. We shall discuss the effects of varying them in Sect. 12.1.41 Specifically, we fiducially adopt td = 3r e f 
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and C ~ 3 for luminous AGNs with peak luminosity L > 10 13 Lq, while td = 0, i.e., the declining phase is almost 
absent for low-luminosity objects. To interpolate continuously between these behaviors, we use a standard erfc-function 
smoothing 


Tp 

T e f 


3 




Aagn \ 

io^W 


(13) 


which is illustrated in Fig. [3] (bottom panel). We note that our results will be insensitive to the detailed shape of the 
smoothing function. The value of £ = 3 is fiducially adopted, since after a time ( td after the peak the BH mass has 
almost saturated to its final value. Results are unaffected by modest variation of this parameter. 

We also fiducially assume that the Eddington ratio Ao of the ascending phase depends on the cosmic time t (or 
redshift z) after 


A oO) = 4 



(14) 


as illustrated in Fig. [3] (top panel). As shown by Lapi et al. (2006, 2014), such moderately super-Eddington values 
at high redshift 2 > 4 are necessary to explain the bright end of the quasar luminosity function (see also Shankar 
et al. 2009, 2013). During the demand-limited, ascending phase of the liglitcurve, Ao exceeds the characteristic 
value Athin ~ 0.3 for the onset of a slim accretion disc (Laor & Netzer 1989). On the other hand, during the 
declining phase of the lightcurve, the Eddington ratio declines almost exponentially, so that after the characteristic 
time Tthin ~ Tp log Ao/Athin the transition to a thin accretion disc takes place. At high redshift where Ao ~ 4, the 
thin-disc regime sets in only after a time Tthin ~ 2.5 rp after the peak, while at low redshift where Ao < 1 it sets 
in about r t hi n ~ 1.2 rp after the peak. We notice that statistically the fraction of slim discs should increase toward 
high- z, as suggested by the data analysis of Netzer & Trakhtenbrot (2014), paving the way for their use as standard 
candles for cosmological studies (Wang et al. 2013c). The time-averaged value of A during the declining phase is, 
to a good approximation, given by (A) ~ Ao (1 — e~^)/( ~ Ao/C, while the time-averaged value during the thin disk 
regime (A) ~ (Athin — Ao e“^)/((/ — log Ao/A t hin) ranges from 0.1 at z < 1 to 0.3 at z > 3. We will see that such values 
(A) averaged over the Eddington distribution associated to the adopted lightcurve reproduce well the observational 
determinations (Vestergaard & Osmer 2009; Kelly & Shen 2013). 

As to the radiative efficiency, we take into account the results of several numerical simulations and analytic works 
(Abramowicz et al. 1988; Mineshige et al. 2000; Watarai et al. 2001; Blandford & Begelman 2004; Li 2012; Begelman 
2012; Madau et al. 2014), that indicate a simple prescription to relate the efficiency e and the Eddington ratio A in 
both slim and thin disc conditions 


Cthin A 
~~2~ e A / 2 - 1 


(15) 


here ethin is the efficiency during the thin disc phase, which may range from 0.057 for a non-rotating to 0.32 for a 
maximally rotating Kerr BH (Thorne 1974). We will adopt ethin = 0.1 as our fiducial value (see Davis & Laor 2011). 
In conditions of mildly super-Eddington accretion with A > a few the radiative efficiency e < 0.3 ethin applies, while 
in sub-Eddington accretion regime with A < 1 it quickly approaches the thin disc value e = ethin- We also take into 
account that along the declining portion of the lightcurve e increases given the almost exponential decrease of A. The 
time averaged values (e) of the efficiency during the declining phase and during the thin disc regime are illustrated in 
Fig. [3] We expect that the redshift dependence of the average efficiency is negligible during the thin disc regime; this is 
in qualitative agreement with the findings by Wu et al. (2013) based on spectral fitting in individual type-1 quasars (see 
also Davis & Laor 2011 for a low -2 determination), and by Cao (2010) based on continuity equation analysis. However, 
we caution the reader that the determination of the radiative efficiency is plagued by several systematic uncertainties 
and selection effects (see discussion by Raimundo et al. 2012). Large samples of AGNs with multiwavelength SEDs 
and BH masses are crucial in fully addressing the issue. 

The final BH mass Mbh is easily linked to the mass at the peak Mbh,p appearing in Eq. m- One has 


rT p+C T n i _ e 

ATbh = / dr '—y Aagn(t') = Mbh.p 

Jo ec 


l + /e^(l-e-<) 

Tef 


(16) 


The correction factor f e takes into account the modest change of the quantity (1 — e)/e along the declining phase. We 
have checked that f e ss 0.8 for any reasonable value of ethin- Notice that at high redshift where Ao ~ 4, the fraction of 
BH mass accumulated in thin-disc conditions is only 5% of the total, while it can be as large as 20% at low redshift 
where Ao ~ 1. This is relevant since most of the BH mass estimates at high -2 are based on single-epoch method, that 
probes the UV/optical bright phase. 

The evolution with the internal time r of the AGN luminosity, mass, and Eddington ratio are sketched in Fig. [2] 
We also schematically indicate with colors the stages (according to the framework described below Eq. fl2l) when the 
galaxy is detectable as a FIR-bright source, and the nucleus is detectable as an X-ray AGN and as an optical quasar. 


2.1.3. BHs: solution 
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Given the lightcurve in Eq. m, the fraction of the time spent by the BH per luminosity bin reads 


E 


d Tj 

dLAGN 


T e f + Tp 

Lagn 


Qh [Lagn < Aagn(ATbh)] 


(17) 


where Lagn (Mbh) is the maximum luminosity corresponding to a final BH mass Mbh> that can be written as 


Lagn(M B h) = 


A Mbhc 

^Edd 


,2 T 


1 + U — (1 - e" c ) 

“Tef 


1 -1 


(18) 


the expression stresses the relevance of the mass accretion during the AGN descending phase. This implies that the 
time spent in a luminosity bin is longer by a factor td than on assuming a simple growing exponential curve, and that 
Eq. (11811 is implicit since td /r e f is itself a function of the luminosity. 

Using Eq. ed in the continuity equation (neglecting dry merging, i.e., no source term) yields 


POO 

Lagn A7(Lagn, t) = / dM B H 9 t 7V(M B H, i) [T e f + r D ] 

-/Mbh(^AGn) 


(19) 


where the minimum final mass that have shone at Lagn is given by the inverse of Eq. m- We proceed by differentiating 
both sides with respect to Lagn and rearranging terms to find 


Lagn <9lagn [L A gn A^Lagn, t)] 

/bH.Lagn r ef + td 


[d t N(M B n,t) Mbh]|m bh (Lagn) ! 


( 20 ) 


in deriving this equation we have defined /bh.l A gn — d log Mbh/^ log Lagn, which is not equal to unity since TD/r e f 
in Eq. m depends on Lagn- 

Finally, we integrate over cosmic time and pass to logarithmic bins. The outcome reads 


N(log 




d t' C^ln Lagn [-N(lpg Lagn )] 


^ -T- , 


( 21 ) 


Note that in practice we have started the integration at z- ln = 10 assuming that the BH mass function at that time 
was negligibly small. This solution constitutes a novel result. In the case when td = 0, and when A and e are constant 
with redshift and luminosity, the above equation reduces to the form considered by Marconi et al. (2004). 


2.1.4. BHs: results 

In Fig. [4] we illustrate our results on the supermassive BH mass function at different representative redshifts. The 
outcomes of the continuity equation can be fitted by the functional shape of Eq. © with Lagn replaced by Mbh , and 
with the parameter values reported in Table 1; the resulting fits are accurate to within 5% in the redshift range from 
0 to 6 and over the BH mass range Mbh from a few 10 7 to a few 10 9 Mg. 

We also illustrate two determinations of the local mass function. One is from the collection of estimates by Shankar 
et al. (2009), that have been built by combining the stellar mass or velocity dispersion functions with the corresponding 
Mbh ^ M* (Haring & Rix 2004) or Mbh — cr (Tremaine et al. 2002) relations of elliptical galaxies and classical bulges. 
The other is the determination by Shankar et al. (2012) corrected to take into account the different relations followed 
by pseudobulges. In addition, we present the determination at z = 0 by Vika et al. (2009) based on an object-by-object 
analysis and on the Mbh — L (McLure & Dunlop 2003) relationship. 

The BH mass function at z ~ 0 from the continuity equation provides and almost perfect rendition of the local 
estimates by Shankar et al. (2009) and Vika et al. (2009) when e t hin = 0.1 is adopted. At 0 w 1 we find a BH mass 
function which is very similar to the local determination. Our result is in good shape with, though on the high side 
of, the determination by Li et al. (2011), based on luminosity (or stellar) mass functions and mild evolution of the 
Mbh — L (or Mbh — M*) relationship. The same also holds at z ~ 2, which is not plotted for clarity. 

At 0 ps 3 we find a BH mass function which at the knee is a factor about 10 below the local data. We are in good 
agreement with the determination by Ueda et al. (2014) based on continuity equation models. This is expected since 
we adopt similar bolometric luminosity functions, and around 0 ss 3 we have similar values of the Eddington ratio and 
radiative efficiency. At z ~ 6 we find a BH mass function which is about 3 orders of magnitude smaller than the local 
data. We compare our result with the estimate by Willott et al. (2010b) in the range Mbh ~ 10 8 — 3 x 10 9 M 0 . This 
has been derived by combining the Eddington ratio distribution from single-epoch BH mass estimates to the optical 
quasar luminosity functions corrected for obscured objects. At the knee of the mass function we find a good agreement 
with our result based on the continuity equation, while at lower masses we predict a slightly higher number of objects. 

The reasonable agreement with previous determinations in the redshift range z ~ 0 — 6 validates our prescriptions 
for the lightcurves, the redshift evolution of Ao(z), and the e — A relation of Eq. ED- Besides, we recall that these were 
already independently tested against the observed fractions of AGNs hosted in sub-mm galaxies and related statistics 
(Lapi et al. 2014; see Sect.[lD. 

Note that during the slim-accretion regime, where most of the BH mass is accumulated, the effective efficiency 
amounts to e < 0.05 given our assumed value ethin ~ 0.1 in Eq. ED. see also Fig. [3] This requires a bit of discussion. 
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In principle, during a coherent disk accretion, the BH is expected to spin up very rapidly, and correspondingly the 
efficiency is expected to attain values e > 0.15 (Madau et al. 2014), corresponding to e t hin ~ 0.3 after Eq. (fl5l) . 
However, such a high value of the efficiency would produce a local BH mass function in strong disagreement with the 
data. This can be understood just basing on the standard Soltan argument. In fact, the BH mass density inferred 
from the AGN luminosity density would amount to pbh ~ 2 x 10 4 (1 — e)/eM 0 Mpc -3 < 10 5 M 0 Mpc -3 . Plainly, 
the z = 0 result would fall short of the local observational determinations, that yields a fiducial mass density of 
Pbh ~ 4.5 x lO 5 Af 0 Mpc -3 (using the Shankar et al. 2009 local mass function). The discrepancy is even worse if 
one considers the local mass function obtained by combining the velocity dispersion or stellar mass function with the 
recently revised Mbh — a or Mbh — M* relations by McConnell & Ma (2013) and Kormendy & Ho (2013), which 
feature a higher overall normalization. 

In App. B we have also tested the relevance of dry merging processes (contributing via the source term in the 
continuity equation) in shaping the BH mass function. At z > 1 BH merging effects are found to be statistically 
negligible (see also Shankar et al. 2009), although smaller mass BHs may undergo substantial merging activity with 
possible impact on the seed distribution (for a review, see Volonteri 2010). At z < 1 our tests indicate that the mass 
function is mildly affected only for Mbh ;£ 10 9 M 0 . 

Thus an average slim-disc efficiency e < 0.05 is required. During the slim-disc accretion, such a low efficiency can be 
maintained by, e.g., chaotic accretion, efficient extraction of angular momentum by jets, or similar processes keeping 
the BH spin to low levels (King & Pringle 2006; Wang et al. 2009; Cao 2010; Li et al. 2012; Barausse 2012; Sesana et 
al. 2014). We also remark that an efficiency e < 0.05 eases the formation of supermassive BHs at very high redshift 

^ 6 , so alleviating any requirement on initial massive seeds (Volonteri 2010). On the other hand, the supermassive 
BH mass function only poorly constrains the values of the BH spin during the final thin-disc phase, which the current 
estimates suggest to be rather high (Reynolds et al. 2013). 

Bolometric corrections and obscured accretion can concur to alleviate the requirement of a low slim-disc efficiency. 
Bolometric corrections are based on studies of SEDs for large samples of AGNs (e.g., Elvis et al. 1994; Richards et 
al. 2006; Hopkins et al. 2007; Lusso et al. 2010, 2012; Vasudevan et al. 2010; Hao et al. 2014). In fact, the SEDs 
depend on the main selection of the objects (e.g., X-ray, UV, optical, IR), possibly on the Eddington ratio (Vasudevan 
et al. 2010; Lusso et al. 2012), and on bolometric luminosity (Hopkins et al. 2007). The recent analysis of Hao et al. 
(2014) finds no significant dependencies on redshift, bolometric luminosity, BH mass and Eddington ratio of the mean 
SEDs for a sample of about 400 X-ray selected type-1 AGNs, although a large dispersion is signalled. A large fraction 
of objects with accretion obscured at wavelengths ranging from X-ray to optical bands has been often claimed, also 
in connection with their contribution to the X-ray background (Comastri et al. 1995). The fraction compatible with 
it at substantial X-ray energies has been recently discussed by Ueda et al. (2014) and properly inserted in our AGN 
bolometric luminosity functions. 

Concerning the overall evolution of the BH mass function, we find that most of the BH mass growth occurs at higher 
redshifts for the more massive objects (see the inset of Fig. HD- The overall BH mass density at z = 0 amounts to 
Pbh ~ 4.5 x 10 5 M 0 Mpc -3 , in excellent agreement with observational determinations. 

In Fig. [5] we show how our results on the mass function depend on various assumptions. The top and middle panels 
illustrate the effect of changing the parameters of the lightcurve: radiative efficiency e, Eddington ratio A, timescale of 
the descending phase 7 d, and duration of the descending phase £. For clarity we plot results only at 2 = 0 and 2 = 3. 
We illustrate our fiducial model, and compare it with the outcome for values of the parameters decreased or increased 
relative to the reference ones. 

To understand the various dependencies, it is useful to assume a simple, piecewise powerlaw shape of the luminosity 
function in the form TV(logLAGN) oc Lagn> with V 1 a t the faint and 77 > 1 at the bright end. Then it is easily seen 
from Eq. m that the resulting mass function behaves as 


A (log M B h) oc 



[1 + (rp/ref) (1 - e c )p 

1 + T D /r ef 


P-^BH 


( 22 ) 


Thus the BH mass function features an almost inverse dependence on e at given BH mass. The dependence on A is 
inverse at the high-mass end, which is mostly contributed by high luminosities where 77 > 1. On the other hand, it is 
direct at the low-mass end, mainly associated to faint sources with 77 < 1. The opposite applies to the dependence on 
To/Vef, since roughly A(logMBH) oc (TD/r e f) 7?_1 . Finally, the dependence on £ is mild, and practically irrelevant for 
£ > 3 since the exponential e~^ in Eq. (1161) tends rapidly to zero. Differences in the results are more evident in the 
2 = 0 than in the 2 = 3 mass function, since this is an integrated quantity, as expressed by Eq. m 

In the bottom left panel we illustrate the effect of changing the optical/X-ray bolometric corrections: the black lines 
refer to our reference one by Hopkins et al. (2007), while the blue and red lines refer to the ones proposed by Marconi 
et al. (2004) and by Lusso et al. (2012), respectively. It is easily seen that the impact on the BH mass function is 
limited, actually well within the uncertainties associated to the input luminosity functions, and to the observational 
determinations of the local BH mass function. 

In the bottom right panel, we illustrate the effect of changing the functional form used to analytically render the 
AGN luminosity functions: the black lines refer to our fiducial rendition via a modified Schechter function (cf. Eq. 0, 
while the green lines refer to a standard double powerlaw representation (e.g., Ueda et al. 2014; Aird et al. 2015). It is 
seen that our results on the BH mass function are marginally affected; this is because both shapes render comparably 
well the input AGN luminosity functions. 
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In Fig. [6] we illustrate the Eddington ratio distribution P(log A|Mbh, z ) associated to the overall adopted lightcurve 
at different redsliift and for different final BH masses. Typically at given redshift and BH mass, the distribution 
features a Gaussian peak at high values of A, and then a powerlaw increase toward lower values of A before an abrupt 
cutoff. The peak reflects the value of A in the ascending part of the lightcurve. Actually since A(r) is constant there, 
the peak should be a Dirac S— function. However, small variations around the central value and observational errors 
will broaden the peak to a narrow Gaussian as plotted here; a dispersion of 0.3 dex has been safely adopted. The 
powerlaw behavior reflects the decrease of A(r) during the declining part of the lightcurve at late times, and the cutoff 
in the distribution mirrors that of the lightcurve. The relative contribution of the Gaussian peak at high A and of the 
powerlaw increase at low A depends on the relative duration of the declining and ascending phases. Thus at given 
redshift, small mass BHs feature a much more prominent peak and a less prominent powerlaw increase relative to high 
mass ones. This is because in small mass objects the descending phase is shorter. At given BH mass, the distributions 
shift to the left, i.e., toward smaller values of A, as the redshift decreases. This is because the initial value Ao(z) 
decreases with redshift, as prescribed by Eq. (El- 

Such a distribution has been computed under the assumption that the overall lightcurve can be sampled. However, 
from an observational perspective, the Eddington ratio distribution is usually determined via single-epoch BH mass 
estimates of type-1 AGNs. This implies that only a portion of the descending phase can be sampled. To ease 
the comparison with observations, we present in the middle and bottom panels of Fig. [G] the expected distribution 
considering only the descending phase (including both the final portion of the slim-disc phase and the while thin-disc 
phase, with A > 0.3), and only the thin-disc phase (i.e., the portion with A < 0.3). The resulting distributions feature 
a powerlaw shape, whose slope depends on the portion of the declining phase that can be effectively sampled: the 
shorter this portion, the steeper the powerlaw. The result is roughly consistent with the observational determinations 
by, e.g., Kelly & Shen (2013), although a direct comparison is difficult due to observational selection effects. In fact, 
different observations are likely to sample diversely the initial part of the declining phase, and this will possibly make 
the expected and the observed distributions even more similar. Note that especially at 2 < 1, BH reactivations, which 
are not included in our treatment (both in terms of lightcurve descriptions and of stochasticity of the events), can 
contribute to broaden the Eddington ratio distribution toward very low values of A < 1CU 2 as estimated in the local 
Universe (e.g., Kauffmann & Heckman 2009; Brandt & Alexander 2015). 

In Fig. [7] we present the AGN duty cycle (<5agn) averaged over the Eddington ratio distribution associated to the 
adopted lightcurve. Specifically, this has been computed as 


{8agn)(Mbh ,t) 


AAGNpOg M B H, t) 

N (log Af B H 5 1) 


dlo g AP(logA|AW)x 


X N(\og TaGN 1 1) Ragn(Mbh.A) 


(23) 


where Tagn(A?bh, A) is given by Eq. (fl8l) . In our approach based on the continuity equation, the duty cycle is 
a quantity derived from the luminosity and mass functions. It provides an estimate for the fraction of active BHs 
relative to the total. At given redshift, the average duty cycle increases with the BH mass, since more massive BHs are 
typically produced by more luminous objects, that feature the descending phase of the lightcurve. On the contrary, 
small mass BHs are originated mainly by low-luminosity objects for which the descending phase is absent. At given 
BH mass, the duty cycle increases with the redshift, essentially because to attain the same final mass, BHs stay active 
for a larger fraction of the shorter cosmic time. This is especially true for BHs with high masses, up to the point that 
they are always active (5agn ~ 1) for z > 3. This agrees with the inferences from the strong clustering observed for 
high-redshift quasars (Shen et al. 2009; Shankar et al. 2010; Willott et al. 2010b; Allevato et al. 2014); we will further 
discuss the issue in Sect. 13.1.11 The increase of the duty cycle with BH mass is consistent with the active fraction 
measured by Bundy et al. (2008), Xue et al. (2010), although the issue is still controversial and strongly dependent 
on obscuration-corrections (see Schulze et al. 2015). On the other hand, we again remark that our approach does not 
include AGN reactivations, which may strongly enhance the duty cycle for low-luminosity objects especially at 2 < 1, 
accounting for the estimates by, e.g., Ho et al. (1997), Green & Ho (2007), Goulding & Alexander (2009), Schulze & 
Wisotzki (2010). 

In Fig. [5] (top panel) we present the AGN Eddington ratio (A) averaged over the lightcurve, computed as 


(X)(M BK ,t) = jv( logMBH ) J dl °S^ ^P(^OgX\M BBl z) N(logL AG N)\L AaN (M BH ,\) ■ (24) 

At given final BH mass, the Eddington ratio decreases with the redshift, as a consequence of the dependence adopted 
in Eq. El- The average values are consistent with those observed for a sample of quasars by Vestergaard & Osmer 
(2009). Note that to take into account the observational selection criteria, we have used the Eddington ratio distribution 
associated to the descending phase, presented in the middle panel of Fig. [ 6 ] 

In Fig. [ 8 ] (bottom left panel) we show the Eddington ratio function, that has been computed as 


N(\ogX,z) 


dlogM B H P(log A|Mbh, ^) AAGN(logM B H, 2) ; 


(25) 


the outcome is mildly sensitive to the lower integration limit, and a value M BB ~ 10 8 Mq has been adopted to compare 
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with data (see Schulze et al. 2015). For the sake of completeness, we present the results when using the Eddington 
ratio distribution associated to the thin disk phase (cf. bottom panel of Fig. [ 6 ]) or to the whole descending phase (cf. 
middle panel of Fig. ED, being the outcome for the overall lightcurve very similar to this latter case. Our results from 
the continuity equation are confronted with the estimates from Schulze & Wisotzki (2010) at z ~ 0, and from Schulze 
et al. (2015) and Nobuta et al. (2012) at 2 ~ 1 — 2, finding a nice agreement within the observational uncertainties 
and the clear systematic differences between datasets. 

In Fig. [5] (bottom right panel) we present a related statistics, i.e., the fraction E(log A|AA*) of galaxies with given 
stellar mass hosting an AGN (active fraction) with a given Eddington ratio. This has been computed simply on dividing 
the quantity P(log A|AAbh, z) NxGN^og Mbh, z) by the number of galaxies N(\ogM+,z) with given stellar mass M* 
(the stellar mass function, cf. Sect. 12.21) . Plainly, AAbh ~ 2 x 10 - 3 AA* must be set to the BH mass corresponding to 
AA*, the result being rather insensitive to the AAbh/AA* ratio adopted; we further take into account a scatter of 0.3 
dex in this relationship, whose effect is to make the active fraction F(logA|AA*) to depend on AA* more weakly than 
the Eddington ratio distribution P(log A|AAbh, z) depends on AAbh- We illustrate the outcome for a range of stellar 
masses from AA* ~ 10 10 5 M 0 to AA* ~ lO 115 Af 0 ; it turns out to be only mildly dependent on AA*, and especially so 
at low redshift z < 2 , as also indicated by current observations. 

In fact, our results can be compared with the observational estimates at z ~ 0 — 2 by Aird et al. (2012) and Bongiorno 
et al. (2012). The latter authors actually provide the active fraction as a function of the observable quantity Lx/M+; 
this can be converted into an Eddington ratio by assuming an X-ray bolometric correction and a value for the AAbh/AA* 
ratio. Bongiorno et al. (2012) suggest an overall conversion factor A ss 0.2 Ax/AA* (here cgs units are used for the 
quantities on the r.h.s.). We also plot their data points when using a conversion A ss 0.08 Lx/M* (corresponding, e.g., 
to a larger ratio AAbh/AA* or a lower bolometric correction), giving more consistency with the determination by Aird 
et al. ( 2012 ). 

All in all, our results from the continuity equation are found to be in good agreement with the observational estimates, 
reproducing their mild dependence on stellar mass and their shape for z < 2 down to an Eddington ratio A ~ a few 
10 -2 . On the other hand, AGNs at 2 < 1 with tiny accretion rates corresponding to A < 10 ~ 2 are likely triggered 
by reactivations, which are not included in our lightcurve, and can contribute to maintain a powerlaw shape of the 
Eddington ratio function and of the active fraction down to A ~ 10 -4 as observed by Aird et al. (2012). 

2 . 2 . The stellar mass function 

Now we turn to the evolution of the stellar mass function from the SFR-luminosity function. 


2.2.1. Stars: SFR function from UV and FIR luminosity 

SFR can be inferred by lines (mainly Lya and Ha) and by continuum emission (mainly UV, FIR, radio and X-ray) 
of star forming galaxies (for a review, see Kennicutt & Evans 2012). The SFR is directly proportional to the UV 
(chiefly far-UV, FUV) luminosity, which traces the integrated emission by young, massive stars. On the other hand, 
even a small amount of dust causes large extinction of the UV emission. The absorbed luminosity is re-emitted at 
longer wavelengths, mostly in the 4 — 1000 /jm interval. Therefore ideal estimates would be based on both the UV 
(L uv) and the FIR (Air) observed luminosities for large galaxy samples at relevant redsliifts. This would allow to 
derive the total luminosity proportional to the SFR 


Asfr = Auv = -^uv + / Air ; (26) 

here the fraction / is meant to subtract from the budget the FIR luminosity contributed by diffuse dust (cirrus) 
absorbing the light from less massive older stars. 

Actually, the SFR functions are inferred from UV or FIR-selected samples. In both cases calibrations and corrections 
come in (Calzetti et al. 2000; Hao et al. 2011; Murphy et al. 2011; Kennicutt & Evans 2012). The calibration constants 
between SFR and luminosity in UV and FIR are practically the same, as expected on energy conservation arguments 
(Kennicutt 1998; Kennicutt & Evans 2012); for FIR luminosity we have 


log 


SFR 

Mq yr - 1 


-9.81 +log/ 


Air 

Am 


(27) 


while for extinction-corrected UV luminosity we have 


log T7r -=T ~ ~ 7A2 ~ °' 4 M %v ~ —9.76 + log uv r , (28) 

Mq yr L,q 

z/uv being the frequency corresponding to 1550 A. @. 

The FIR luminosity ascribable to the diffuse dust emission (cirrus) depends on several aspects such as stellar content 
(mass, age and chemical composition), dust content and spatial distribution. The cirrus emission is characterized 
by dust temperature lower than the emission associated with star formation in molecular clouds (Silva et al. 1998; 
Rowlands et al. 2014). There are local galaxies with quite low SFR, whose FIR luminosity is dominated by cirrus 
emission. E.g., Hao et al. (2011) found 1 — / ~ 0.5 for a sample of nearby starfornhng galaxies with SFR AA* < 30 Af 0 


8 Some UV data are given at a restframe wavelength A different from 1550 A; Eq. 11281 still holds provided that on the right hand side 
the correction — log A/1550 A is added. E.g., for A = 1350 A the correction amounts to 0.06 and the zero point calibration becomes —7.36. 
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yr . However, the fraction 1 — / strongly reduces with increasing star formation (e.g., Clemens et al. 2013). For 
strong local starbursting galaxies with M* > 100M 0 yr -1 and Lir > 10 12 L 0 , such as Arp 220, we get 1 — / < a few 
percent (Silva et al. 1998; Rowlands et al. 2014). Hereafter we will assume that / = 1 for Lir > 10 12 Lq and that at 
such large luminosities the SFR can be soundly inferred from the FIR luminosity. 

At low luminosity, the SFR is better estimated from UV emission. For this purpose it is essential to allow for dust 
absorption, that may drastically reduce the UV luminosity to a few percent or even less of its intrinsic value. When 
only UV data are available, the correlation between the UV slope /3 and the IRX ratio Tir/Tuv is largely used to 
infer the dust attenuation (Meurer et al. 1999). While initially proposed only for low redshift galaxies, the method 
has been tested and applied also to high redshifts (Reddy et al. 2010; Overzier et al. 2011; Hao et al. 2011; Bouwens 
et al. 2013, 2015). However, for large values of the slope /? and of the attenuation, the spread around the correlation 
becomes huge (Overzier et al. 2011; Reddy et al. 2010) and the estimate of attenuation becomes quite uncertain even 
in local samples (e.g., Howell et al. 2010). On the other hand, the estimate of attenuation for UV-selected samples 
is less dispersed for galaxies with SFRs M* < 1 Mq yr -1 . In such instances the correction to UV luminosity is more 
secure and relatively small on average (Bouwens et al. 2013). In fact, this is also suggested by the UV attenuation 
inferred by combining the Ha attenuation and the Calzetti extinction curve (Hopkins et al. 2001; Mancuso et al. 
2015). 

Given all these considerations, we build up the overall SFR-luminosity Usfr function as follows. We start from 
the luminosity function at different redshifts observed in the FIR band by Magnelli et al. (2013), Gruppioni et 
al. (2013), Lapi et al. (2011), and in the UV band by Bouwens et al. (2015), Oesch et al. (2010), Reddy & 
Steidel (2009), and Wyder et al. (2005). The data are reported in Fig. [9] In passing note that the SFR and 
the SFR-luminosity Lsfr scales on the upper and lower axis have been related assuming the approximate relation 
logSFR/M 0 yr" 1 a: —9.8 + logisFR/Ug, and so the number density per unit SFR or per unit luminosity is the same. 

For the FIR samples we assume / = 1, while for the UV samples at redshift z > 2 we have exploited the dust 
correction suggested by Meurer et al. (1999) and Bouwens et al. (2013, 2015). At z < 2 the attenuation has been 
kept to the values found by Bouwens et al. (2013) for z « 2.5 galaxies. This assumption at z < 1 produces an UV 
attenuation somewhat in between those proposed by Wyder et al. (2005) and Cucciati et al. (2012), and the one 
proposed by Hopkins et al. (2001). However, we stress that for galaxies with Lyv 10 10 Lq the correction is anyway 
smaller than a factor ~ 2. 

Fig. [9] shows that at any redshift we lack a robust determination of the SFR-luminosity function at intermediate 
luminosities; this occurs for two reasons: first, UV data almost disappear above Luv ~ 10 11 Lq (see also Reddy 
et al. 2010) because of dust extinction, while FIR data progressively disappear below L\jy « 10 12 Lq because of 
current observational limits. Second, the UV correction for Tuv ^ 10 10 Lq or intrinsic SFR M* > 1 Mq yr -1 becomes 
progressively uncertain, as discussed above. 

To fill in the gap, we render the overall SFR distribution with a continuous function, whose shape is basically 
determined at the bright end by the FIR data and at the faint end by the UV data. Specifically, we exploit the same 
modified-Schechter functional shape of Eq. ©, with Tagn replaced by Tsfr and the parameter values reported in 
Table 1. The UV data at the faint and FIR data at the bright are smoothly connected by our analytic renditions at 
various redshifts. We also illustrate an estimate of the associated ler uncertainty. In the inset we plot the evolution 
with redshift of the SFR-luminosity density and the contribution to the total by specific luminosity ranges. 

It happens that our rendition of the data closely follows the model proposed by Mao et al. (2007) and Cai et al. 
(2014), wherein the extinction is strongly differential with increasing SFR (and gas metallicity). In such models, the 
faint end of the UV luminosity function at high redshift is dictated by the rate of halo formation, while the bright end 
is modeled by the dust content in rapidly starforming galaxies. At z > 6 reliable statistics concern only UV-selected 
galaxies endowed with low SFR. At high luminosity we have extrapolated the behavior from lower redshift z < 4, 
finding a good agreement with the model proposed by Cai et al. (2014). This extrapolation implies, at z > 6, a 
significant fraction of dusty galaxies with SFR M* > 10 2 Mq yr” 1 , which are missed by UV selection. Clues of such 
a population are scanty, but not totally missing. Riechers et al. (2014) detected a dust obscured galaxy at £ ~ 6.34 
with SFR M* ss 2900 Mq yr" 1 , and Finkelstein et al. (2014) a second one at z ~ 7.51 with SFR M* ss 300 M 0 yr” 1 . 
The large SFR end at z > 6 will be probed in the near future by ALMA and JWST observations. 

In passing, we have also reported the extrapolation of the SFR-luminosity function to z = 10 (cyan line in Fig. [9]). 
It is interesting to compare this with the recent estimates from UV observations by Bouwens et al. (2015). At L\jy « 
10 9 7 Lq the extrapolation matches the observed number density around 10 3 Mpc 3 . For smaller Tuv ~ 10. 9 7 Lq we 
remark that the slope of the luminosity function is highly uncertain; data extrapolation suggests a slope in the range 
from —1.65 to —2, as illustrated by the cyan shaded area. At the other end, for Luv ~ lO lo4 Z/ 0 , the extrapolated 
number density is around 10" 4 Mpc" 3 , a factor around 3 times larger than that observed in the UV. This possibly 
suggests that dust already at z ~ 10 affects the UV data toward the bright end, as it happens at lower redshift. 

2.2.2. Stars: lightcurve 

There are three time-honored assumptions regarding the behavior of the SFR as a function of galactic age: expo¬ 
nentially increasing (up to a ceiling value), exponential decreasing, constant. 

Here we specialize to a very simple, constant lightcurve, motivated by the recent FIR data from the Herschel 
satellite concerning high-redshift, luminous starbursting galaxies, and their physical interpretation on the basis of the 
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BH-galaxy coevolution model by Lapi et al. (2014). Hence we adopt 

£SFR(T|M*, t) = ft* M* T < Tb urs t 


— 0 T /> Tburst 


(29) 


where k* is a dimensional constant converting SFR into bolometric luminosity. For a Chabrier IMF we have k* ft 
2.5 x 10 43 yr erg s _1 /M 0 ft! 6.5 x 10 9 yr Lq/Mq (see Sect. I2.2.T1) . The constant SFR M* = AAburst/Tburst represents 
an average over the fiducial period of the burst Tburst, with the total mass of formed stars amounting to A/^burst- 

Since the more massive stars restitute most of their mass to the ISM, the total amount of surviving mass is Af* = 
(1 — 1Z) M*,burst, where 7Z is the restituted fraction that depends on the IMF and on the time elapsed from the burst. 
For the Chabrier IMF the mass in old, less massive stars approaches to 1 — 1Z ft 0.5 when the time elapsed is larger 
than a few Gyrs. Since we shall exploit the continuity equation also at relatively short cosmic times at z > 1 we adopt 
the value 1 — 1Z ft 0.6 corresponding to Tburst ~ 1 Gyr (see below). 

The most recent observations by ALMA have undoubtedly confirmed that the SFR in massive high-redshift galaxies 
must have proceeded over a timescale around < 0.5 Gyr at very high rates > a few xlO 2 Af 0 yr” 1 under heavily 
dust-enshrouded conditions (e.g., Scoville et al. 2014, 2015, their Table 1). The observed fraction of FIR detected host 
galaxies in X-ray (e.g., Page et al. 2012; Mullaney et al. 2012b; Rosario et al. 2012) and optically selected (e.g., Mor 
et al. 2012; Wang et al. 2013b; Willott et al. 2015) AGNs points toward a SFR abruptly shutting off after this period 
of time. In the analysis by Lapi et al. (2014) this rapid quenching is interpreted as due to the energy feedback from the 
supermassive BH growing at the center of the starbursting galaxy. In the first stages of galaxy evolution the BH is still 
rather small and the nuclear luminosity is much less than that associated to the star formation in the host. The SFR 
is then regulated by feedback from SN explosions, and stays roughly constant with time, while the AGN luminosity 
increases exponentially. After a period < 1 Gyr in massive galaxies the nuclear luminosity becomes dominant, blowing 
away most of the gas and dust from the ambient medium and hence quenching abruptly the star formation in the host. 
On the other hand, longstanding data on stellar population and chemical abundances of galaxies with final stellar 
masses Af* < 10 10 Af 0 indicate that star formation have proceeded for longer times regulated by supernova feedback 
(see reviews by Renzini et al. 2006; Conroy 2013; Courteau et al. 2014; and references therein). 

On this basis, we adopt a timescale for the duration of the starburst given by 


Tburst (t) = 1 Gyr 



-3/2 


1 + 2 erfc 



Asfr 
10 105 L 



(30) 


the dependence on the cosmic time mirrors that of the dynamical/condensation time, in turn reflecting the increase 
of the average density in the ambient medium. In addition, the erfc-smoothing connects continuously the behavior 
for bright and faint objects expected from the discussion above. We tested that our results are insensitive to the 
detailed shape of the smoothing function. At high redshift, as noted by Lapi et al. (2014), such a timescale is around 
15 — 20 e—folding time of the hosted BH (i.e., < 0.5 — 1 Gyr). The luminosity scale 3 x 1O 1O L 0 corresponds to SFR 
Af* ft 5 Af 0 yr” 1 . 


2.2.3. Stars: solution 

Given the lightcurve in Eq. (BSD. the fraction of the time spent per luminosity bin reads 


E 


d Ti 

d+SFR 


Tburst 5d [Lsfr — Asfr(A1*)] 


(31) 


with Lsfr(M*) = k* M^burst/Tburst the SFR-luminosity associated to the final stellar mass Af*; the Dirac delta- 
function Sd(-) specifies that, since the lightcurve is just a constant, the luminosity associated to a stellar mass Af* 
must be in the luminosity bin d+sFR- 

Using this expression in the continuity equation Eq. m without source term yields 


' burst J sfr 

where the final stellar mass that have shone at Lsfr is given by 

Af*(L S FR, t) = a-^SFR Tburst _ 


(32) 


(33) 


In deriving these equation, we have used dAf*/dLsFR = f*,L SFR A/*/Lsfr- On the same line of Sect. l2.1.3l we integrate 
over cosmic time, and turn to logarithmic bins. The outcome reads 


_/V(log Af*,t) 



Af(logLsFR) 

f*,L sfr T'burst 




(34) 
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This solution constitutes a novel result. Note that our approach exploits in the continuity equation the full SFR- 
luminosity function, and is almost insensitive to initial conditions; in these respects it differs from the technique 
developed by Leja et al. (2015; see also Peng et al. 2010) to evolve the stellar mass function backwards from z < 2 
basing on the observed SFR—M* relationship and the starforming fraction. 

Interestingly, if Tb urs t is independent of TsfRi a Soltan-type argument can be extended to the stellar content. It can 
be easily found on multiplying Eq. (13411 by M* and integrating over it, to obtain 

f dM*M* AT(M*,i) = - f d t' f dLsFR^SFR AT(Lsfr,0 ; (35) 

i/ (/ 0 »/ 

comparing with the classic expression for the BH population, it is seen that the role of the efficiency combination 
(1 — e)/ec 2 ~ 7 x 10 -14 (1 — e)/e yr _1 Mq/Lq is played by the quantity (1 — 7 £)/k* ~ 9 x 10 -11 yr _1 Mq/Lq, which 
mainly depends on the IMF (here the constant refers to the Cliabrier’s IMF). 

In passing, we notice that for conventional IMFs most of the stellar mass in galaxies resides in stars with mass 
< 1 Mg. Since these stars emit most of their luminosity in the near IR, the galaxy stellar mass M* can be inferred by 
the near-IR luminosity functions. At variance with the BH case, the so called ’remnants’ are not dark but luminous red 
stars. This provides a significant vantage point to estimate the mass function of these ’remnants’. In fact, the stellar 
mass function N(M+) is worked out via the statistics of the stellar luminosity function N(L+), not to be confused with 
the SFR-luminosity function N(Lg FR ) used above. 

2.2.4. Stars: results 

In Fig. [TO] we illustrate our results on the stellar mass function at different representative redshifts. The outcomes 
of the continuity equation can be fitted with the functional shape of Eq. |9]) with Lagn replaced by M*, and the 
parameter values given in Table 1. The resulting fits are accurate within 5% in the redshift range from 0 to 6 , and in 
the stellar mass range M* from a few 10 9 to a few 10 12 Mq. 

The outcome at z ~ 0 is compared with the determination of the local mass function by Bernardi et al. (2013). The 
outcomes at z « 1 and z « 3 are compared with the determinations by Santini et al. (2012) and Ilbert et al. (2013), 
while the result at z ~ 6 is compared with the measurements by Stark et al. (2009) and Duncan et al. (2014). The 
results of the continuity equation and the observational estimates at different redshifts are in very good agreement, 
considering the associated uncertainties and systematic differences among different datasets. 

Concerning the overall evolution, the high-mass end of the mass function is mainly built up at z > 1.5, while 
the low-mass end is still forming at low z. The inset shows the progressive buildup of the stellar mass density as a 
function of redshift. The global stellar mass density at z = 0 reads p* « 3 x 10 8 Mq Mpc -3 , in good agreement with 
observational determinations, and a factor about 10 3 above the total BH mass density. The stellar mass densities at 
z ~ 1 is already very close to the local value. 

In Fig. HH we show how our results on the stellar mass function depends on the parameters of the lightcurve: the 
timescale of burst duration Tb urs t and the adopted IMF. To understand the various dependencies, it is useful to assume 
a simple, piecewise powerlaw shape of the luminosity function in the form N(\ogLs FF ) oc Lgp R , with 77 < 1 at the 
faint and 77 > 1 at the bright end. Then it is easily seen from Eq. (l34l) that the resulting stellar mass function behaves 
as 

AC (log M*) oc T bnrltVM- ri . (36) 

Thus the stellar mass function features an almost direct dependence on Tburst at the high-mass end, which is mostly 
contributed by high luminosities where 77 > 1. On the other hand, the dependence is inverse at the low-mass end, 
mainly associated to faint sources with 77 < 1. The dependence on the IMF is related to the ratio (1 — 7£)//c*; e.g., 
passing from the Chabrier to the Salpeter (1955) IMF, the ratio increases by a factor of 2. More significant variations 
are originated when passing from Chabrier to a top-heavy IMF (e.g., Lacey et al. 2010) which implies the ratio to be 
reduced by a factor ~ 8 . 

We have also tried to parameterize the stellar lightcurve with a decreasing or increasing exponential like Tsfr oc 
e -T / T *; the solution of the continuity equation in these instances can be derived on the same route used for BHs. The 
net result is that to reproduce the observed stellar mass function at different redshifts the typical timescale of the 
exponential r* must be of the order of Tb urs tj he., the lightcurve is required to be approximately constant over such a 
timescale as we have indeed assumed. 

Fig. fl2l shows the average duty cycle (<5sfr) of star formation in galaxies. In analogy to Eq. |0, this has been 
computed as (^sfr)(T/*, z) = N [logM*(Lg FR ), z ]/N (log M+,z), where the relation between Tsfr and M* is given by 
Eq. (1291) . At z > 1 the duty cycle is almost unity, reflecting the build up of the stellar mass function in real time. On 
the other hand, at z < 1 the duty cycle progressively drops down, dramatically for stellar masses M* > a few 10 11 Mq. 
This is because the mass added by in situ star formation becomes negligible. 

Two related outcomes presented in the Appendices are extremely relevant in this context, the first concerning dust 
formation, and the second concerning the role of dry merging. In App. C we highlight the fundamental role of the 
dust, by confronting our fiducial result with that derived basing on the UV-selected luminosity functions. Fig. IC2I 
directly shows that the UV-selected luminosity function, even corrected for dust extinction, produces a stellar mass 
function much lower than the observed one for M* >2x 10 10 Mq at any redshift z < 6 . In particular, we stress that 
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our extrapolated FIR portion of the SFR-luminosity function at z ~ 6 is validated by the good comparison with the 
stellar mass function observed around that redsliift. This implies that massive galaxies formed most of their stars in 
a dusty environment. We expect a large fraction of massive galaxies to be already passively-evolving (i.e., with quite 
low SFR and ‘red’ colors) at z > 1, as indeed increasingly observed even at substantial redsliift z ~ 3 (Man et al. 
2014; Marchesini et al. 2014). 

The point is strengthened by Fig. 1131 which shows an estimate of (actually an upper bound to) the dust ‘formation’ 
time Tdustj computed multiplying the starformation timescale iburst by the ratio of the UV-selected to the total SFR- 
luminosity functions. At z ~ 6, galaxies with SFRs M* s=s 100 M 0 yr -1 and final stellar masses M* > 3 x 10 10 M 0 have 
a dust formation time of Tdust ~ 3 x 10 7 yr, implying a quite rapid metal/dust enrichment. Interestingly, this is much 
shorter than the fiducial time ss 15r e f «a few 10 8 yr to grow the hosted final BH mass (cf. Eq. [12]). Therefore, most 
of the BH growth must occur in dusty galaxies (e.g., Mortlock et al. 2011). At redsliift z > 2 — 3 the constraints on 
t dust for strongly starforming objects stays almost constant. In moving toward lower redshift z < 2 the dust formation 
time becomes shorter, even for moderately starforming objects with SFR M* < 30 M 0 yr -1 . This can be interpreted 
as star formation episodes mainly occurring within dust-rich molecular clouds, within galaxies already evolved as to 
the chemical composition of their ISM. 

In App. B we investigate the impact of dry merging on the evolution of the stellar mass function. In this context it 
is worth stressin g tha t the effect of dry merging is negligible at redshift z > 1 and it can play some role only at lower 
redshift (see Fig. IB2I) . These outcomes statistically ascertain that most of the stellar content in massive galaxies with 
M„ >3x 10 10 M 0 is formed in situ. However, we caution that the observed stellar mass function cannot currently be 
assessed for M* > 3 x 10 11 M 0 given the substantial systematic uncertainties in the data (see discussion by Bernard! et 
al. 2013), and that the role of dry mergers can be of some relevance in the growth of such extremely massive galaxies 
(see Liu et al. 2015; Shankar et al. 2015). 

All in all, we stress the capability of the continuity equation in reconstructing the star formation history in the 
Universe from the past SFR activity. 


3. ABUNDANCE MATCHING 

Having obtained a comprehensive view of the bolometric luminosity and mass functions for stars and supermassive 
BHs at different redshift, we now aim at establishing a link among them and the gravitationally dominant dark matter 
component. To this purpose, we exploit the abundance matching technique, a standard way of deriving a monotonic 
relationship between galaxy and halo properties by matching the corresponding number densities (Vale & Ostriker 
2004, Shankar et al. 2006, Moster et al. 2010, 2013; Behroozi et al. 2013). 

When dealing with stellar or BH mass M, we derive the relation M(Mu. z) with the halo mass Mh by solving the 
equation (e.g., White et al. 2008; Shankar et al. 2010) 


/log M 


dlogM' IV (log M', z) = 




( \/2 (7 log 


M 


(37) 


which holds when a lognormal distribution of M at given Mh with dispersion eri og m is adopted. In the above expression 
we have defined <7i og m = cr\ og M /’/x with /x = d log M/d log Mr. On the basis of the investigation by Lapi et al. (2006, 
2011, 2014) on the high-redshift galaxy and AGN luminosity function, we expect the Mbh — Mh correlation to feature 
a quite large scatter <Tio g M B H ~ 0-4 dex, while a smaller value eri og M* ~ 0.15 dex is expected for the correlation with 
the stellar component. We shall compare the correlations M — Mh obtained when such values for the scatter are 
considered with those obtained by assuming one-to-one relationships, i.e., taking <ri og m = 0. 

In Eq. J37D the quantity JV/logMn) is the galaxy halo mass function (GHMF), i.e., the mass function of halos hosting 
one individual galaxy. We do not simply rely on the overall halo mass function (HMF), because we aim at obtaining 
relationships valid for one single galaxy, not for a galaxy system like a group or a cluster. In a nutshell, we build up 
the GHMF on correcting the overall HMF from cosmological ./V—body simulations, by adding to it the contribution 
of subhalos, but by probabilistically removing from it the contribution of halos corresponding to galaxy systems. We 
defer the reader to App. A for the detailed description of this procedure. The resulting GHMF is plotted at different 
redsliifts in Fig. [[4j we stress that the determination of the GHMF as a function of redshift constitutes in its stand a 
novel result. The outcomes can be fitted with the functional shape of Eq. ([9]) with Lagn replaced by Mh, and with 
the parameter values reported in Table 1. The resulting fits are accurate within 5% in the redshift range from 0 to 10. 

In the same Figure we also compare the GHMF to the overall HMF. The difference between the two, i.e. the galaxy 
system halo mass function at z = 0, is compared with local data to cross-check our approach. At the bright end the 
GHMF drastically falls off, so that even at 2 < 1 galactic halo masses of ss 10 14 M 0 are very rare, since these masses 
pertain to galaxy systems. These findings agree with galaxy-galaxy weak lensing measurements (Kochanek et al. 2003; 
Mandelbaum et al. 2006, van Uitert et al. 2011, Leauathaud et al. 2012, and Velander et al. 2014) and dynamical 
observations in nearby galaxies (Gerhard et al. 2001, Andreon et al. 2014; see also the review by Corteau et al. 2014). 

The same abundance matching technique may also be applied to the stellar or AGN bolometric luminosity L, looking 
for a relation L(Mh, 2 ) specifying the typical luminosity to be expected in a halo of mass Mh at given redshift 2 . 
However, when dealing with luminosities, one has to take into account that galaxies and AGNs shine only for a fraction 
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of the cosmic time. In practice, we use a modified abundance matching of the form 


d log L 


'log L 


, N(log L',z) 
(5) x t 


+ OO 


1 


d log d+N(log Mg, z) - erfc 


f log[M H (A)/M^] 

<Tlog L 


(38) 


where ( S) x t is the duty cycle 6 averaged over the lightcurve multiplied by the cosmic time t, and N(log Mu, z) is 

the formation rate of galactic halos computed according to Lapi et al. (2013). 

3.1. Abundance matching results 

We turn to present the results of the abundance matching technique among various statistical properties of BH, 
galaxies, and host halos. Analytic fits to such outcomes can be found in App. D and Table 2. 


3.1.1. BH vs. halo properties 

In Fig.[15]we show the relationship between the final BH mass Mbh and halo mass Mg from the abundance matching 
technique, at different redshifts. 

Since we are comparing BH and halo statistics at fixed z, the resulting relationship constitutes a snapshot, that can 
be operationally exploited in numerical works to properly populate halos at the reference redshift. On the other hand, 
the evolution of BHs and halos due to accretion is expected to modify, though on different timescales, the relation as 
the cosmological time passes. For example, if the cosmological growth of halos is dominant, then the relation would 
shift along the Mh axis. The relationship at a subsequent redshift takes into account such an evolution, although the 
number of evolved BHs and halos is generally subdominant with respect to the newly formed objects. 

The top panel of Fig. IT5l shows the results when a one-to-one (i.e., no scatter) Mbh — Mg relationship is assumed, 
while bottom panel shows the resulting average relationship when a Gaussian distribution in Mbh at given Mg with a 
scatter of 0.4 dex is adopted. The presence of scatter is particularly relevant at high redshift. Assuming a one-to-one 
relationship would yield at z « 6 average BH masses Mbh > 10 10 Mq within halos of Mh >5x 10 12 Mg, much larger 
than at z ~ 3. This would imply a significant change in the physical mechanisms establishing the Mbh — Mh relation 
over a relatively short timescale of ~ 1 Gyr. In the presence of scatter instead such large BH masses constitute only 
extreme instances, relative to much smaller average values Mbh ~ 10 9 M 0 , not very different from the lower redshift 
ones. In this scenario such peculiar instances are precisely those picked by current observations at high-redshift, which 
are biased toward the more luminous AGNs powered by the more massive BHs. Thus the one-to-one relationship 
offers a view of the observed properties at the given redshift, while the average relationship (with scatter) is helpful to 
provide a physical interpretation. With scatter included, taking into account the considerable uncertainties, one can 
estimate that the logarithmic slope of the average relationship at z > 1 is around Mbh oc M ^ , so encompassing the 

range suggested for AGN feedback processes (Silk & Rees 1998; Fabian 1999; King 2005; for a recent review see King 
2014). The average relationship is practically unchanged within the uncertainties over the range z ~ 1 — 6. Plainly, at 
z = 0 we find very good agreement with the relation inferred from the BH mass function by Shankar et al. (2009). 

In Fig. [16] we show the relationship between the AGN luminosity Tagn and halo mass Mh, both with and without 
scatter. Concerning the scatter, the same comments of the previous Figure apply. The flattening in the relation toward 
lower redshift is mainly driven by the evolution of the AGN luminosity function, especially at the bright end. The 
one-to-one relationship, together with the duty cycle, is required to properly populate halos and derive the clustering 
properties of AGNs. 

In Fig. |TT] we show the luminosity- and BH mass-averaged AGN bias as a function of redshift z. This has been 
computed as follows: we start from the linear halo bias model 6 (Mh,z) including excursion set peak prescriptions 
(Lapi & Danese 2014) for halos of mass Mh at redshift z (see also Sheth et al. 2001). Then we associate to each halo 
mass Mh an AGN of luminosity Lagn as prescribed according to the Lagn — Mh one-to-one relationship discussed 
above. Finally, we compute the luminosity-weighted bias as a function of redshift 

L, , _ IZ L min dlo & L N (log L AGN ,z)b(L AGN, z) ^ 

KZ) Q Li n in cnogL AGN N(logL AGN ,z) ’ ( j 

where L m - m is a minimum bolometric luminosity. The same procedure can be followed to obtain the AfsH-averaged 
bias through the one-to-one Mbh — Mh relation and the average over the BH mass function. 

The resulting bias as a function of redshift compares well with the observational data points from large optical and 
X-ray survey samples. Typically, the optical data refer to quasars with luminosities Lagn a few 10 12 Lq, while X-ray 
data refer to AGNs with Lagn ^ a few 10 11 L 0 ; however, the selection of the datasets reported in the plot are diverse, 
and the reader is deferred to the original papers for details. Note that while at z > 2 X-ray selected AGNs appear 
to be less clustered than optical quasars, the opposite holds true at low z < 2. This fact is somewhat puzzling since 
X-ray AGNs feature generally lower bolometric luminosities, and is often interpreted in terms of a different accretion 
mode becoming dominant at low z (e.g., sporadic reactivation episodes in place of continuous accretion, see discussion 
by Allevato et al. 2011, 2014). As a reference, the halo bias 6 (Mh, z) for various Mh is also shown. It is evident that 
typical host halos feature Mh > 10 12 M 0 with a clear tendency for more massive halos to host more luminous AGNs 
and more massive BHs. In the inset it is seen that even the mild trend of the bias with luminosity at given redshift 
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2 ~ 2 from optical surveys is reproduced. On the other hand, the dependence on luminosity is expected to significantly 
increase at higher z > 4. 

We stress that the clustering properties constitute a byproduct of our approach, and the comparison with observations 
validate our results on BH mass function and duty cycle (see also Shankar et al. 2010). Note that past studies (Martini 
& Weinberg 2001) have instead exploited the clustering properties to constrain the AGN duty cycle. In comparing 
with previous works related to the AGN bias (e.g., White et al. 2008; Hopkins et al. 2007; Wyithe & Loeb 2009; 
Bonoli et al. 2010; Shankar et al. 2010), a few remarks are in order: (i) we stress that our adoption of the GHMF 
in place of the routinely-used HMF appreciably improves the agreement with observations of the bias for luminous 
AGNs/massive BHs at z > 3; (ii) we confirm that values A > a few at 2 > 3, implying a quite rapid growth of the BH 
during the ascending portion of the AGN lightcurve, are also required to meet the observational constraints; (iii) we 
find that the weak dependence of the bias on luminosity at z ~ 2 is rather insensitive to the presence of the descending 
portion of the AGN lightcurve, that we recall is instead indicated in luminous objects by the observed fraction of 
starforming hosts in optically-selected quasars (see Sect. 12.13 . 

3.1.2. Stellar vs. halo properties 

In Fig. HH]we show the relationship between the final stellar mass M„ vs. the halo mass Mh, for different redshift. 
The result at 2 = 0 is compared with the relationship inferred from the local stellar mass function by Bernardi et 
al. (2013). We find a good agreement within the associated uncertainties. The M* vs. Mh at given redshift can 
be described by a powerlaw with slope around 1 at the high-mass end, then steepening for halo masses Mh < a few 
10 12 M 0 . The presence of the scatter around 0.15 dex does not affect appreciably the correlation. 

At 2 > 1 the statistics of both stellar and halo masses are dominated by newly-created objects, so that the evolution 
in both masses of the older individuals is irrelevant. From this perspective, the little if no evolution of the M* — Mh 
relationship can be interpreted in the light of similar, in-situ processes regulating the star formation at different 
redsliifts (Moster et al. 2013). This may be seen more clearly in the inset, showing the efficiency M*//bMn for the 
conversion into stellar component of the original baryon content within the halo /b Mh, having adopted a cosmic initial 
baryon to DM ratio /b = 0.2. The efficiency rises from values < 10% for halo masses Mh < 10 11 M 0 to a roughly 
constant values < 25% around halo masses Mh ~ a few 1O 12 M 0 . All in all, the star formation process in halos is 
highly inefficient. From a physical point of view, this is usually interpreted in terms of competition between cooling 
and heating processes. In low-mass halos, the latter is provided by energy feedback from SN explosions. In massive 
halos, cooling rates are not significantly depressed by SN feedback, and the star formation can proceed at much higher 
levels until the AGN attains enough power to quench it abruptly. 

At 2 < 1 the interpretation is more complex, since the statistics of stars and halos are no longer dominated by 
newly-formed objects and significant evolution in one of the two components may occur. Specifically, for high masses 
the halo evolution dominates, and the M* — Mh evolves shifting toward higher halo masses at almost constant stellar 
mass; contrariwise, for small masses, stellar mass evolution dominates over the halo’s, and the relationship shifts 
upward at almost constant halo mass. 

In Fig. [T9] we present a comparison of our M* — Mh relationship at z = 0 with literature determinations. Our result 
when the GHMF is exploited for the abundance matching (same as in previous Figure) can be directly compared with 
the determination by Shankar et al. (2006) based on the same abundance matching technique. The difference is mainly 
due to the dynamical adopted by Shankar et al. in building the stellar mass function from the If—band 

luminosity function. 

On the other hand, our result when the overall HMF is adopted can be directly compared to the determinations 
based on the abundance matching by Behroozi et al. (2013) and by Moster et al. (2013). These are quite similar to 
ours at the low-mass end, while appreciably steeper at the high-mass end (see also Kravtsov et al. 2014; Shankar et 
al. 2014), where the Bernardi et al. (2013) stellar mass function we adopt contains relatively more objects. 

These results based on the overall HMF can also be directly compared with the data from gravitational lensing 
measurements in groups and clusters of galaxies by Han et al. (2014), Velander et al. (2014), and Mandelbaum et al. 
(2006). The agreement is very nice. Note that since gravitational lensing probes the mass projected along the line of 
sight, it is sensitive to the presence of groups and/or clusters surrounding the individual galactic halo. 

In Fig. [20] we show the relationship between the luminosity Lsfr associated to the SFR vs. the halo mass Mh, 
for different redshifts. The presence of the scatter around 0.15 dex only marginally affects the average relationship. 
We show both the outcome based on the overall SFR-luminosity function, and the one based on the dust-corrected 
UV luminosity function only. This has been determined by matching the GHMF with the SFR-luminosity function of 
UV-selected galaxies corrected for dust extinction (see App. C and Fig IC1I) . It is evident that the typical UV data 
substantially underestimate the luminosities associated to the SFR. We stress once again that FIR data are crucial for 
an unbiased view of the star formation process. 

In Fig.[2n]we also plot the relationship expected at 2 « 10, although we caution that for halo masses Mh < 10 11 M 0 
the relationship strongly depends on the faint-end slope of the luminosity function. To illustrate the variance, we plot 
as a lower bound the relation corresponding to a faint-end slope —1.65, and an upper bound corresponding to —2 
(Bouwens et al. 2015; see also Sect. I2.27TI) . The latter instance is required to keep the Universe ionized out to 2 < 8.8, 
corresponding to an electron-scattering optical depth r es ~ 0.066 as recently measured by the Planck Collaboration 
(2015). Our SFR vs. Mh relationship suggests that this can be afforded by galaxies starforming at rates > 10~ 2 M 0 
yr -1 , with UV magnitudes Muv ^ —13 hosted within halos of masses Mh > 1O 9 M 0 (see also Cai et al. 2014). 

In Fig. [21] we show the galaxy bias, both luminosity (or SFR)-averaged and stellar mass-averaged, for different 
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values of minimum SFR or M*. These quantities have been computed following the same procedure for the AGN bias 
as described in Sect. 13.1.11 For reference we also report the halo bias for different halo masses. It is seen that the 
bias computed from the abundance matching reproduces very well the determination at different redshifts for various 
populations of objects. In particular, UV-selected objects like Lyman Break Galaxies and Lyman-a emitters feature 
low stellar masses M* < 1O 9 M 0 and SFRs less than a few M 0 yr -1 , while FIR-selected objects are associated to 
much more violent SFRs > 10 2 M 0 yr -1 and constitute the progenitors of massive galaxies with final stellar content 
M* > 10 11 Mq. 


3.1.3. SFR and sSFR vs. stellar mass and redshift 

In Fig.[55]we plot the cosmic specific (sSFR) defined as the ratio between the SFR density psfr = / d log M+ M* AT(log M+) 
and the stellar mass density p* = f d log M* M* N (log M*). The resulting cosmic sSFR reflects the behavior for typical 
SFR and stellar masses at the knee of the corresponding distributions, and it includes all galaxies, even the passively 
evolving ones (see also Madau & Dickinson 2014). 

We report both the outcome based on the overall SFR-luminosity function, and the one based on the dust-corrected 
UV luminosity function only. It is apparent that the latter case underestimates the cosmic sSFR at any redshift (cf. 
Wilkins et al. 2008). We also illustrate the result by Madau & Dickinson (2014), which is similar to ours up to 2 ~ 2, 
and then approaches the UV-inferred result. As a matter of fact, their cosmic star-formation history at z > 3 is based 
on UV data (see their Fig. 9). 

The reported observational estimates refer to galaxy samples selected with different criteria. Specifically, at z > 3 
they mainly refer to UV-selected samples. This explains why they are better reproduced by our results for the UV 
dust-corrected case. On the other hand, at z < 1.5 they are mainly based on UV+near-IR data with ongoing star 
formation inferred from 24 /jm or radio fluxes. In this redshift range the sSFR estimated from the ratio psfr/p* 
lies below most of the data points, because it includes an increasing fraction of objects in passive evolution, while 
observations refer to starforming galaxies only (see discussion by Madau & Dickinson 2014). On the other hand, Ilbert 
et al. (2015) report values of the sSFR closer to the ratio psfr/p +, but a factor of ~ 2 — 3 lower than previous estimates 
in literature. The authors attribute this difference to their more accurate treatment of the selection effects, leading to 
inclusion of galaxies with lower sSFR, and to their more accurate statistics. 

In Fig. [231 we show the relationships between the SFR and the stellar mass M*, at different redshifts; this is often 
referred to as the ‘main sequence’ of starforming galaxies (e.g., Elbaz et al. 2011; Rodighiero et al. 2011). Note that 
the outcome is obtained by matching the abundances of two observable quantities like the SFR-luminosity and stellar 
mass functions (the halo mass is bypassed), including the star formation duty cycle. From this point of view, the 
outcome is only mildly dependent on assumptions on the star formation lightcurve and timescales. As in the previous 
Figure, we report the outcome from the abundance matching based on the overall SFR-luminosity functions, and the 
one based on dust-corrected UV luminosity functions only. 

We compare the abundance matching result with the recent observational estimates by Rodighiero et al. (2014), 
Speagle et al. (2014), Salmon et al. (2015), Renzini & Peng (2015) based on large samples of individual measurements 
of SFRs and stellar masses. We stress that, especially at z < 1, determinations of the main sequence by various 
authors differ, mainly because of the way galaxies are selected as being starforming (see discussion by Renzini & Peng 
2015). Further observations and analysis are needed to fully assess the main sequence, both regarding the the overall 
normalization and the slopes and the high and low mass end (that can even be different, see Whitaker et al. 2014). 

At 2 > 1 our results from the abundance matching based on FIR+UV luminosity functions well agree with the 
estimates by Rodighiero et al. (2014) and Speagle et al. (2014) based on multiwavelength observations of galaxy 
samples. At 2 ss 6 our result from the abundance matching based on UV luminosity function is in excellent agreement 
with the data for UV-selected galaxies by Salmon et al. (2015). 

On the other hand, for z < 1 our results appear to be at variance with the observational determinations for 
stellar masses M* < 5 x 1O 1O M 0 . However, in this range the results from the abundance matching becomes loosely 
constraining, because of the large uncertainties introduced by the flatness of the stellar mass function (cf. Fig. fTIil). 

This suggests that the stellar mass and SFR luminosity function may not sample the same galaxy population at their 
respective faint end. Nevertheless, at high masses the abundance matching technique is consistent with current data, 
and actually extends the main sequence in a range where determinations from individual measurements are still scanty. 

3.1.4. BFl mass vs. stellar mass 

In Fig.[24]we illustrate the relationship between the BH mass and stellar mass at different redshifts. The computation 
is performed by the abundance matching of the BH and stellar mass function from the continuity equation, thus 
bypassing the halo mass. We show results both for the one-to-one case (top panel), and when a Gaussian scatter of 0.3 
dex between Mbh and M* is assumed (bottom panel). The presence of the scatter is increasingly relevant at higher 
redshift in biasing observations toward extreme values of the Mbfi/M* ratio (shown in the inset). It is worth noticing 
that the evolution of the relationship and hence of the Mbr/M+ ratio is quite small for z < 3, at variance with the 
claims by some authors (Peng 2007; Jahnke & Maccio 2011). This signals once again that the BH and stellar mass 
growth occurs in parallel by in-situ accretion and star formation processes. 

Our results at 2 = 0 agree with the relations inferred from the abundance matching of the local determinations for 
the stellar and BH mass functions by Bernardi et al. (2013) and Shankar et al. (2009). Our findings are in very good 
agreement with the individual determinations of BH and stellar masses based on dynamical measurements by Haring 
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& Rix (2004). On the other hand, Kormendy & Ho (2013) propose a relation which is systematically higher by a factor 
ft 2.5. The Soltan argument would then imply an extremely high final BH mass density, and in turn a value e < 0.02 
of the average efficiency during the slim-disc regime (see Sect. H). 

In Fig. [25] we illustrate the evolution with redshift of the mass density in DM halos, stars, and BHs. The stellar 
mass density closely mirrors that of galactic DM halos, because the star to DM mass ratio (i.e., the star formation 
efficiency) stays roughly constant with redshift for typical galaxies at the knee of the mass function (see Fig. fl8l) . 
On the other hand, for z < 2 the stellar mass density progressively differs from that of the overall halo population 
(including galaxy groups/clusters). Once again this strengthens the point that star-formation at high redshift occurs 
via in-situ processes within the central regions of galactic halos. The stellar mass density is a factor about ~ 30 — 50 
lower than the galactic halo mass density, reflecting the inefficiency of galaxy formation due to feedback processes, as 
discussed in Sect. 13.1.21 

The BH mass density has a considerably different shape, that toward higher z progressively steepens relative to the 
galactic halo and stellar mass density. This is due to two effects: (i) the number density of halos able to host massive 
BHs declines rapidly; (ii) the time needed to grow massive BHs becomes comparable with the age of the Universe, 
so making apparent the delay of about a few 10 8 yrs between the BH and stellar formation. In the inset we show 
that the observed density ratio between the SFR and the AGN luminosity attains a minimum around z ~ 1.5 and it 
stays almost constant toward lower redshift. This is because both luminosity densities decline in parallel (cf. insets in 
Figs. [T|and [9]). The same trend also applies for the corresponding mass density ratio. 

4. SUMMARY 

We have investigated the coevolution of galaxies and hosted supermassive black holes throughout the history of the 
Universe by a statistical approach based on the continuity equation and the abundance matching technique. Our main 
results are the following: 

• We have demonstrated that the local supermassive BH mass function and the stellar mass functions at different 
redshift can be reconstructed from the SFR and AGN luminosity functions via a continuity equation approach 
without source term. This implies that the buildup of stars and BHs in galaxies occurs mainly via local, in-situ 
processes, with dry mergers playing a marginal role at least for stellar masses M* <3x 10 11 M 0 and BH masses 
A/bh < 10 9 ATq, where the statistical data are more secure and less biased by systematic errors. 

• As to the AGN/BH component, our analysis nicely reproduces the observed Eddington ratio function and the 
observed fraction of galaxies with given stellar mass hosting an AGN with given Eddington ratio (see Sect. 12.1 Hi) . 
Such an agreement strongly suggests that the fraction of AGNs observed in slim-disc regime increases with 
redshift, and that anyway most of the BH mass is accreted in such conditions. 

• We have inferred relationships between the stellar, BH, and DM components of galaxies at various redshifts. 
These imply that stellar and AGN feedback cooperate with gas cooling in the star formation process within halos, 
whose binding energy at formation is the most relevant feature. Specifically, in low-mass halos SN explosions keep 
star formation low on long timescales, while in massive halos star formation can proceed at much higher levels until 
the AGN quenches it abruptly. These relationships between galaxy/BH and halo properties constitute testbeds 
for galaxy formation and evolution models, and can be operationally implemented in numerical simulations to 
populate dark matter halos or to gauge subgrid physical prescriptions. Duty cycles for both the AGN and the 
stellar components are derived, and found to be close to unity at high-redshift. 

• We have derived the bias as a function of redshift and luminosity, both for the AGN and for various galaxy 
populations. The clustering properties constitute a byproduct of our approach, and the nice agreement with 
observations validate our results on BH and stellar mass functions, and related duty cycles from the continuity 
equation. 

• The specific SFR increases with redshift at least up to z ~ 6. In the range z > 1 the results from the abundance 
matching technique agree with the so called ‘main sequence’ of starforming galaxies, although we underline that 
the comparison with observations critically depend on sample selection. For z < 1 the results from abundance 
matching are reliable for stellar masses M* > 5 x 10 10 Mq, where they are consistent and actually extend the 
observational determinations in a range where individual measurements are still scanty. 

• We show how strongly the presence of the dust affects the view of the star formation process in galaxies with 
SFRs M* > 10 Mq yr” 1 at any redshift, even the quite large ones. In fact, we have shown that dust is formed on 
a timescale which is only a small fraction of the burst duration. Such a behavior is also mirrored in the estimated 
cosmic SFR and sSFR density. 

• The low efficiency < 20% in star formation elucidates that a fraction > 50%, up to ~ 70% depending on mass, 
of the gas associated to a galaxy halo is always in warm/hot form. 

• The BH to stellar mass ratio evolves mildly at least up to z < 3, signaling that the BH and stellar mass growth 
occurs in parallel by in-situ accretion and star formation processes. 
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The marginal role of dry merging and the inefficiency of star formation imply that galaxy formation is basically a 
process inherent to the inner regions of halos, where most of the gas mass resides. 

These evidences strongly add motivation to the the development of hydrodynanrical simulations at very high spatial 
resolution, which allow detailed studies of small-scale gravitational instabilities connected to gas cooling and condensa¬ 
tion, star formation, BH accretion, and associated feedback processes (e.g., Ceverino et al. 2015; for a comprehensive 
review Bournaud 2015). Our main results are listed in Table 3, where we also recall their location in the paper, and 
cross-reference to the corresponding sections and figures. 

From the technical point of view, the novel achievements of the present work can be summarized as follows: 

• We have presented an analytical solution of the continuity equation for BHs that holds under quite general as¬ 
sumptions, including a redshift/luminosity dependence of the Eddington ratio, radiative efficiency, and lightcurve 
timescales. 

• We have developed the continuity equation for the stellar component, solving it under quite general assumptions 
about the lightcurve shape and timescales. 

• We have provided a continuous rendition of the overall SFR function, interpolating between the UV data at the 
faint and the FIR data at the bright end. A posteriori, our approach is validated by the agreement of the stellar 
mass function via continuity equation with the observational determinations over the redsliift range 2 ~ 0 — 6. 

• We have developed a procedure to derive the galaxy halo mass function at different redshifts. This can be 
implemented in halo occupation distribution models. 

• We have generalized the abundance matching technique to deal with relationships between luminosity and mass, 
by considering the duty cycle of BHs and star formation in galaxies. 

We stress that the added value of continuity equation and abundance matching is to provide largely model- 
independent outcomes, which must be complied by detailed physical models. 

Finally, two remarks are in order. As to the AGN/BH component, large samples of AGN with multiwavelength 
SEDs are crucial in testing the statistics of the slim-disc fraction and in measuring the associated radiative efficiency 
(cf. Raimundo et al. 2012). As to the stellar component, our analysis allows to extrapolate the SFR, stellar mass, and 
sSFR functions to higher redsliift, yet unexplored but within the reach of future instrumentations like ALMA , JWST 
and SKA. In particular, a crucial point will be to estimate the bright end of the SFR-luminosity function at 2 > 4, to 
obtain direct constraints on the timescale of dust formation in high-redsliift galaxies. 
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APPENDIX 

APPENDIX A: GALACTIC HALO MASS FUNCTION 

In this Appendix we detail our procedure to derive the galactic halo mass function, i.e., the mass function associated 
to halos hosting one individual galaxy. The computation actually includes two steps: (i) we account for the possibility 
that a halo contains various subhalos; (ii) we probabilistically exclude halos corresponding to galaxy systems rather 
than to individual galaxies. 

Our starting point is the subhalo mass function, as recently determined by Jiang & van den Bosch (2014). The 
distribution of subhalos with mass between in and in + dm in a halo of mass Mjj at redsliift 2 can be well fitted by 
the function 

N(logif>) = 7 ^“ e -/3 ^ In 10 , (Al) 

where if) = m/M h- Actually if m is taken as the subhalo mass at the infall time, the resulting unevolved subhalo mass 
function is universal for any mass Mh and as such described by the parameter set [ 7 , a, /3, w] = [0.22, —0.91,6.00,3.00]. 
This is plotted in Fig. Al. 

However, we are more interested in taking to as the mass of the surviving, self-bound entity at redshift z, which is 
reduced with respect to that at accretion due to mass stripping and dynamical friction. The resulting evolved subhalo 
mass function is then described by the same functional shape in Eq. (Al) but with modified parameter set [ 7 , a, /3, w] = 
[0.31 f s , —0.82, 50.00,4.00] which depends on the host halo mass and redshift through the quantity f s . The latter may be 
determined as follows: first one obtains the half-mass redshift 20.5 solving ^(20.5) = < 5 C ( 2 ) + 1.19 \Z<j 2 (M-r/ 2 ) — ct 2 (M h ), 
6 c (z) being the linear threshold for collapse at redshift 2 , and a 2 (M) the mass variance at the scale M. Then one 
computes N T = J dt/rd yn (t), Tdyn being the halo dynamical time. Finally, f s = 0.3563 ALT 0 6 — 0.075 holds. The 
outcome is illustrated for different redshift and host halo masses in Fig. Al. 
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Now we can compute the overall subhalo contribution to the halo mass function in the mass bin between Mh and 
Mh + dMn as 

pOO 

AUbH(logM H ) = / dlogMn A^H(logMn) A’(logV’)|v.=M H /M^ , (A2) 

Jo 

where AnQogMH) is the standard halo mass function (see Sheth & Tormen 1999; Tinker et al. 2008). Thus the total 
halo + subhalo mass function just reads 

Adr+subH (log Mr) = Ah (log Mh) + A su bH(log Mh) . (A3) 

In Fig. A2 we plot at different redshift the halo mass function, the overall subhalo mass function, and the total halo 
plus subhalo mass function. It is easily seen that the subhalo contribution is almost negligible for any redshift in the 
mass range of interest for this work. 

Now we turn to compute the probability distribution for a given halo to contain one individual galaxy. The first 
step is to obtain the halo occupation number (HON), i.e., the average number of subhalos inside a host halo of mass 
Mh; it writes 

(A)(M h , z) = f dlog-0 N(log ip) . (A4) 

J log m m i n / Mh 

Here m m ; n represents a minimum mass for subhalos, required to avoid the divergence in the above integral. This will 
be set by comparison with numerical simulations and observational datasets. The resulting HON as a function of Mh 
and redshift, for different minimum subhalo masses m m i n is illustrated in Fig. A3. For high Mh m m ; n the HON is 
well fitted by a power-law with logarithmic slope < 1, going into an abrupt cutoff for masses Mh < 3 — 5m m i n - 
The HON represents the average number of subhalos inside a host halo, but we need instead the probability dis¬ 
tribution P(N\(N)) of having N subhalos given the average number (A)(Mh, z). Numerical simulations and HOD 
models aimed at reproducing various galaxy observables (Zeliavi et al. 2005, 2011; Zheng et al. 2007, 2009; Tinker et 
al. 2013) indicate that such a distribution is well approximated by a Poissonian. Then one can easily compute the 
cumulative probability P{< N\(N)) of having less than N subhalos. This reads 


P(<N\(N)) 


t(n + i,(n)) 

N\ 


(A5) 


where T(a, x ) = f°° d tt a 1 e 1 is the incomplete complementary F-function, x is the floor function (the closest integer 
lower than x), and n\ = 1 x 2 x ... x n the factorial function. We stress that in such a probability the dependence on 
host halo mass and redshift are encased into the HON (iV)(M H , z). 

Finally, the galaxy halo mass function can be computed as 

AcHMF(logM h ) = AH+ su bH(logM H ) x P(< N = 1|(A)) . (A6) 

The outcomes at different redshifts, having adopted a minimum satellite mass of m m ; n = 10 11 Mg, are illustrated in 
Fig. [14] of the main text. With respect to the full halo + subhalo mass function, the GHMF features a cutoff for host 
halo masses Mh > 1 — 3 x 1O 13 M 0 , more pronounced at lower redshift. This is because the probability of hosting 
subhalos (hence more than one galaxy) increases strongly for large masses. In other words, such halos are more likely 
to host a galaxy group or cluster than an individual galaxy. 

We stress that both a minimum mass of ~ a few 10 11 M 0 for satellite halos (corresponding to the adopted m m ;„), 
and a maximal value ~ a few 10 13 M 0 for a halo to host an individual galaxy (corresponding to the resulting cutoff 
in the GHMF) are strongly indicated by galaxy weak lensing observations (Mandelbaum et al. 2006, van Uitert et al. 
2011, Leauathaud et al. 2012, and Velander et al. 2014). Furthermore, such maximum galactic halo masses are also 
strongly suggested by dynamical observations of gas and stars in nearby galaxies (see Gerhard et al. 2001, Andreon 
et al. 2014; see also the review by Corteau et al. 2014). 

Finally, to provide a further observational test, we have computed the group and cluster mass function by subtracting 
the GHMF from the overall halo + subhalo mass function. This represents the halo mass function associated only to 
galaxy systems, and as such is expected to show a cutoff for halo masses < 1O 13 M 0 . The result at z = 0 is plotted 
as a dotted line in Fig. [14] of the main text, and compared with the determinations by Boehringer et al. (2014) from 
X-ray observations of groups and clusters and by Martinez et al. (2002) from optical observations of loose groups; the 
agreement is quite impressive. We notice that the behavior in the range of poor clusters and groups is particularly 
sensitive to the parameter m m i n . Thus the agreement with the data is another indication that the adopted value 
around 10 11 M 0 is appropriate. 

The overall halo and the galaxy halo mass functions can be fitted with the functional shape of Eq. ([9]) with Tagn 
replaced by Mh, and the parameter values given in Table 1. The resulting fits are accurate within 5% in the redshift 
range from 0 to 10 and halo mass range from 10 to 14 for the GHMF and from 10 to 16 for the HMF. 


APPENDIX B: DRY MERGERS 

Many recent works (e.g., Shankar et al. 2009, 2013; Moster et al. 2010, 2013) have shown that the role of dry mergers 
(i.e., addition of the whole mass content in stars or BHs of the merging halos without contributing significantly to 
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star formation or BH accretion) in building up the BH/stcllar mass functions is less important than accretion/in-situ 
star formation at z > 1. This is simply because the evolutionary times associated to mergers are much longer than 
those associated to the in-situ BH/stellar mass growth. On the other hand, at z < 1 the situation is expected to 
reverse. This is because the cold accreting or starforming gas within the DM halo gets progressively exhausted or is 
ejected/heated by the energy feedback from supernovae or from the AGN itself. In fact, the continuity equation with 
accretion only yields little evolution of the mass functions from z ~ 1 to z ~ 0 (cf. Figs. HI and flTTl and related insets). 
Thus the low redshift z < 1 evolution could be in principle more affected by dry merging. This is a hotly debated 
issue in the literature, especially in relation to the size evolution of massive, passively-evolving galaxies (e.g., Naab et 
al. 2009; Fan et al. 2010; Nipoti et al. 2012; Shankar et al. 2015; Kulier et al. 2015). 

In this Appendix we highlight the impact of dry mergers on the supermassive BH/stellar mass functions at z < 1. 
We start from the observed mass functions at z ~ 1, then evolve them down to redshift zero by taking into account 
both dry mergers and accretion in the continuity equation. The effect of dry mergers is evaluated numerically with 
a mid-point scheme computation that divides the overall timegrid in sufficiently small steps St, and then evolves the 
mass function at each timestep + according to 


N(log M,U + St) = N(log M,ti) + ^ N(log M/2,U) St - V N(\og M,ti) St , (Bl) 

where V is the probability of dry mergers. We adopt the common simplifying assumption that dry merging of the 
associated stellar and BH components follows halo mergers of given mass ratio. We base on the DM merging rates 
provided by Stewart et al. (2009) via high-resolution A"—body simulations, and write 

st n — n 0 - 72 

n> 0 - 0.02 ±- (l + *) X1 1 ^ , (B2) 

where £ specifies the mass ratio above which mergers are considered; thus V(> 0.5) is the probability of major mergers, 
while V(> 0.1) — V(> 0.5) is that of minor mergers. 

The results on the BH and stellar mass functions are illustrated in Fig. Bl and B2. The impact of dry mergers on 
the mass functions is apparent only at the high mass end. Dry mergers increase moderately the space densities of BHs 
with mass Mbh > 10 9 Mq and boost that of stellar masses M* > 10 12 M 0 , ranges where data are still statistically 
uncertain and/or affected by large systematics. 

Specifically, assuming that a dry merger follows any DM halo merger (either major or minor) yields a local BH 
mass function still consistent with data, even considering the uncertainties on the bolometric corrections in converting 
from luminosity to mass (cf. Sect. [UTTTTl) . On the other hand, while major dry merger produce a stellar mass function 
consistent with the data (see also Liu et al. 2015), this is not the case for minor dry mergers. All that implies that the 
addition of stellar mass by (minor) dry mergers following the DM halos’ must be only partial, possibly depending on 
mass ratio, orbital parameters, tidal stripping, and structural properties (see Naab et al. 2009; Krogager et al. 2014). 

APPENDIX C: THE COMPLEMENTARITY OF UV AND FIR DATA 

In this Appendix we stress the importance of the FIR, in addition to the UV, data in probing the star formation 
process in high-redshift galaxies. 

To this purpose, we present in Fig. ICll the SFR luminosity function estimated on the basis of dust-uncorrected 
UV data, dust-corrected UV data, and UV+FIR data. It is evident that, even when dust-corrected according to the 
prescriptions described in Sect. 12 .2.11 UV data strongly undersample the bright end of the luminosity function. For 
example, at £ ~ 3 the number of sources with M* ~ 300 Mq yr — , which is not an extreme but rather a typical value, 
is estimated to be 10~ 4 Mpc -3 from UV+FIR data, while it is inferred to be < 10~ 6 Mpc -3 from dust-corrected UV 
data, and would be < 10 _1 ° Mpc -3 from dust-uncorrected UV data. We stress that especially at z > 1.5, the dust 
corrections routinely applied to the UV data are unable to fully account for the population of strongly starforming 
galaxies seen in the FIR band. 

In Fig. lC2l we illustrate the stellar mass function obtained via the continuity equation from the above input luminosity 
functions. We keep the same lightcurve of our fiducial model, that for UV-bright, low-luminosity objects yields a 
duration of the burst already close to the Hubble time. As it can be seen, when basing only on UV data (even if 
dust-corrected) the high mass end of the resulting stellar mass function is substantially underpredicted relative to the 
FIR+UV results (which well reproduces observational estimates, see Fig. [TUI) at any redshift. Note that this mismatch 
can hardly be recovered by mass additions from dry merging events, since a factor 10 in mass is needed from z ~ 3 to 
0 . 


APPENDIX D: ANALYTIC FITS TO ABUNDANCE MATCHING RELATIONSHIPS 

Here we provide analytic fits to the relationships derived from the abundance matching technique. To fit a relation 
of the form Y — M we adopt a double powerlaw shape: 


M 


_M b (z)_ 


*(*) 


+ 


M 


lM b (z)\ 


+z)' 


Y(M,z) = N(z) x 


(Dl) 
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with 6 = — 1 for a convex or 9 = +1 for a concave relationship. 

The normalization log_/V( 2 ), the mass of the break log Mb(z), and the characteristic slopes a(z) and u(z) evolve 
with the redshift according to the same parametrization 

p{z) = Po + k p i x + k p2 X 2 + k P 3 X 3 (D2) 

with 

*= log (rvi) (D3) 

and z o = 0.1. The parameter values are reported in Table 2. 

These expressions can be exploited to interpolate and/or extrapolate the relationships all the way from 2 « 0 to z ~ 6. 
Interpolation is helpful to produce mock galaxy and AGN/BH catalogs that can be used to compute gravitational 
lensing effects, to investigate clustering properties, to gauge sub-grid physics in numerical simulations, and to design 
observational setups. On the other hand, extrapolation is particularly helpful to obtain specific predictions in redshift 
and mass ranges not currently probed by the data but within the reach of upcoming experiments. 
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Figure 1. The bolometric AGN luminosity function A(logLAGN) at redshift z = 0 (orange), 1 (red), 3 (green), and 6 (blue). Optical 
data are from Richards et al. (2006; filled circles), Croom et al. (2009; filled squares), Masters et al. (2012; filled triangles), Ross et al. 
(2012; filled stars), Fan et al. (2006; filled pentagons), Jiang et al. (2009; filled reversed triangles), Willott et al. (2010a; filled diamonds); 
X-ray data are from Ueda et al. (2014; open squares), Fiore et al. (2012; open stars), and Aird et al. (2015, open diamonds). The optical 
and X-ray luminosities have been converted to bolometric by using the Hopkins et al. (2007; see their Fig. 1) corrections, while the number 
densities have been corrected for the presence of obscured AGNs accord ing to Ueda et al. (2003, 2014). The solid lines illustrate the 
analytic rendition of the luminosity functions as described in Sect. 12.1.11 and the hatched areas represent the associated uncertainty; the 
cyan line is the extrapolation to z = 10 plotted for illustration. The inset shows the AGN luminosity density as a function of redshift, for 
the overall luminosity range probed by the data (solid line with hatched area), and for AGN bolometric luminosity logT agn/ i n the 
ranges [11,12] (dot-dashed line), [12,13] (dashed line), [13,14] (dotted line). 
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Figure 2. Time evolution of the luminosity (top panel), mass (middle panel), and Eddington ratio/specific SFR (bottom panel) normalized 
to the value at time of the AGN luminosity peak. Solid lines refer to AGN-related, and dashed lines to star-formation related quantities. 
The orange area highlights the stage when the galaxy is starforming and appears as a FIR, bright source (orange); the red and bl ue a reas 
highlight the stages when the BH is detectable as an X-ray AGN and as an optical quasar, respectively. See main text below Eq. 112> for 
details. 
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Figure 3. Top panel: The adopted Eddington ratio (magenta lines) and radiative efficiency (green line) as a function of redshift. The 
values in the ascending, demand-limited phase (solid lines) and the time-averaged values during the descending, supply-driven phase 
(dashed lines) and during the thin-disc regime (dotted lines) are also shown. Bottom panel: the characteristic timescale t^/tq f of the 
AGN descending phase (magenta line) and the duration Tb urst of the stellar burst (green lines) at redshift z = 1 (dashed), 3 (solid), and 6 
(dotted) as a function of the peak AGN and of the SFR luminosity, respectively. 
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Figure 4 . The supermassive BH mass function JV(logMBH) as a function of final BH mass Mbh- Results from the continuity equation 
(see Sect. [2JL31 at redshift z = 0 (orange), z = 1 (red), 3 (green), and 6 (blue) are plotted as solid lines, with the hatched areas illustrating 
the associated uncertainty; the cyan line is the extrapolation to z = 10 plotted for illustration. The dark grey shaded area illustrates the 
collection of estimates by Shankar et al. (2009) built by combining the stellar mass or velocity dispersion function with the Mbh — M* or 
Mbh — & relations of elliptical galaxies; the light shaded area is the determination by Shankar et al. (2012) corrected to take into account 
the different relations followed by pseudobulges. The orange circles illustrate the determination at z = 0 by Vika et al. (2009). The red 
dashed area illustrate the determination at z ~ 1 by Li et al. (2011), the green dashed area shows the range of models by Ueda et al. (2014) 
at z ~ 3, and the blue dashed area the estimate by Willott et al. (2010b) at 2 : ~ 6. The inset shows the BH mass density as a function of 
redshift computed from the continuity equation, for the overall mass range (solid line with hatched area), and for BH masses log Mbh/Mq 
in the ranges [6,7] (dot-dashed line), [7,8] (dashed line), [8,9] (dotted line). The grey shaded areas illustrate the observational constraints 
from the above z = 0 mass function by Shankar et al. (2009, 2012). 
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Figure 5. Comparison plot showing the dependence of the supermassive BH mass function on various assumptions; for clarity only results 
at 2 : = 0 (solid lines) and at z = 3 (dashed lines) are plotted. In the top and middle panels, we show the effects of changing the values of 
the parameters describing the AGN lightcurve. The black lines are for our fiducial values, the red and blue lines refer to values decreased 
or increased of the amount specified in the legend, and the green lines to constant values in redshift and luminosity. In the bottom left 
panel we illustrate the effect of changing the optical/X-ray bolometric corrections: the black lines refer to our reference one by Hopkins et 
al. (2007), while the blue and red lines refer to the ones proposed by Marconi et al. (2004) and by Lusso et al. (2012), respectively. In 
the bottom right panel, we illustrate the effect of changing the functional form used to analytically render the AGN luminosity functions: 
the black lines refer to our fiducial rendition via a modified Schechter function (cf. Eq. [9]) , while the green lines refer to a double powerlaw 
representation. 
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Figure 6. The Eddington ratio distribution P(log A|Mbh ? z) associated to the overall lightcurve (top panel), only to the descending phase 
(middle panel), and only to the thin-disc phase (bottom panel), at different redshift z = 0 (orange), 1 (red), 3 (green), and 6 (blue) and 
for different BH masses Mbh — 10 6 (dotted), 10 7 (dashed), 10 8 (solid), and 10 9 Mq (dot-dashed); for clarity the results relative to masses 
10 6 , 10 7 and 10 9 M© are plotted only at 2 : = 0. 













CONTINUITY EQUATION AND ABUNDANCE MATCHING 


33 


Average duty cycle 



z 


Figure 7. The average AGN duty cycle (<?agn) 38 a function of redshift z , for different BH masses Mbh = 10 6 (dotted), 10 7 (dashed), 
10 8 (solid), and 10® Mg (dot-dashed). The inset illustrates the AGN duty cycle as a function of the BH mass at different redshift z = 0 
(orange), 2 = 1 (red), 3 (green), and 6 (blue). 
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Figure 8. Top panel: the average Eddington ratio (Aagn) as a function of redshift z, for different BH masses Mbh — 10 6 (dotted), 
10 7 (dashed), 10 8 (solid), and 10 9 M@ (dot-dashed); the cyan shaded area covers the range of measured values by Vestergaard Sz Osmer 
(2009). Bottom left panel: the Eddington ratio function at redshift z = 0 (orange), 1 (red), and 2 (green) associated to the thin-disk phase 
(solid lines) or to the descending phase (dashed lines; the outcome when considering the overall lightcurve is very similar); observational 
estimates at z ~ 0 are from Schulze Sz Wisotzki (2010; orange diamonds), at z ~ 1 — 2 from Schulze et al. (2015; red circles for VVDS, 
squares for zCOSMOS, and triangles for SDSS) and from Nobuta et al. (2012; red stars for SXDS). Bottom right panel: the fraction of 
galaxies of given stellar mass (solid lines for M* 10 10-5 Mq and dotted lines for M* ~ 10 11,5 Mq) hosting an AGN with given Eddington 
ratio at redshift z = 0 (orange), 1 (red), and 2 (green); data are from Aird et al. (2012; orange squares) and from Bongiorno et al. (2012; 
red and green circles), where the latter have been converted with the relation A « 0.2 Lx/M* (filled circles) or A ~ 0.08 Lx /M* (empty 
circle), see text for details. 
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Figure 9. The SFR-luminosity function TV(logLsFR) at redshift z = 0 (orange), 1 (red), 3 (green), and 6 (blue), vs. the bolometric 
luminosity Lsfr associated to the SFR (lower axis) and vs. the SFR (upper axis). Infrared data are from Magnelli et al. (2013; filled 
circles), Gruppioni et al. (2013; filled squares), and Lapi et al. (2011; filled stars); UV data (dust corrected, see text) are from Bouwens 
et al. (2015; open circles), Oesch et al. (2010; open squares), Reddy &; Steidel (2009; open stars), Wyder et al. (2005; open diamonds), 
Finkelstein et al. (2014; open triangles), Cucciati et al. (2012; open reversed tria ngles) , Weisz et al. (2014; pentagons). The solid lines 
illustrate the analytic rendition of the luminosity functions as described in Sect. 12.2.11 and the hatched areas represent the associated 
uncertainty; the cyan line is the extrapolation to 2 = 10 plotted for illustration, with the shaded area representing the uncertainty on the 
faint-end slope. The inset shows the SFR-luminosity density as a function of redshift, for the overall luminosity range probed by the data 
(solid line with hatched area), and for bolometric luminosity logLgFR /-^0 i n the ranges [ 10 , 11 ] (dot-dashed line), [ 11 , 12 ] (dashed line), 
[12,13] (dotted line). 
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Figure 10. The stellar mass funct ion AT(logM*) as a function of the (survived) final stellar mass M* in solar units. Results from the 
continuity equation (see Sect. 12.2.31) at redshift z = 0 (orange), z = 1 (red), 3 (green), and 6 (blue) are plotted as solid lines, with the 
hatched areas illustrating the associated uncertainty; the cyan line is the extrapolation to z = 10 plotted for illustration. High redshift data 
are from Ilbert et al. (2013; filled stars), Santini et al. (2012; filled squares), Stark et al. (2009; filled diamonds) and Duncan et al. (2014, 
filled pentagons). Local data at z = 0 are from Bernardi et al. (2013): filled circles with errorbars illustrate their fiducial measurements 
with the associated statistical uncertainty, while the shaded area shows the systematic uncertainty related to light profile fitting. The inset 
shows the stellar mass density as a function of redshift computed from the continuity equation, for the overall mass range (solid line with 
hatched area), and for stellar masses log M^/Mq in the ranges [9,10] (dot-dashed line), [10,11] (dashed line), [11,12] (dotted line). The 
grey shaded area illustrates the observational constraints from the z = 0 mass function. 
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Figure 11. Comparison plot showing the dependence of the stellar mass function on the parameters of the assumed stellar lightcurve; for 
clarity only results at 2 : = 0 (solid) and z = 3 (dashed lines) are plotted. In the top panel, the black line is our fiducial model, while the 
red and blue lines refer to values decreased or increased relative to the reference one; the green lines refers to a constant (in redshift and 
luminosity) value. In the bottom panel, the black line refers to our fiducial Chabrier IMF, while the colored lines are for Kennicutt (1983; 
red), Salpeter (1955; blue) and a top-heavy (Lacey et al. 2010; green) IMF. 
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Figure 12. The average stellar duty cycle (<5sfr) as a function of redshift z, for different stellar masses M* = 10® (dotted), ! 0 1 0 (dashed), 
10 11 (solid), and 10 12 A/q (dot-dashed). The inset illustrates the duty cycle as a function of the stellar mass at different redshift z = 0 
(orange), 2 = 1 (red), 3 (green), and 6 (blue). 
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Figure 13. An estimate of (actually an upper bound to) the dust formation timescale as a function of the SFR-luminosity (lower scale) 
and of the SFR (upper scale) at redshifts z = 1 (red), 3 (green) and 6 (blue), computed from dust-corrected UV data (dot-dashed lines), 
and dust-uncorrected UV data (dotted lines); for comparison the timescale of the burst duration is also shown (solid lines). 
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Figure 14. The galaxy halo mass function N(log Mr) (solid lines) at redshift z = 0 (orange), 1 (red), 3 (green), 6 (blue) and 10 (cyan), 
vs. the halo mass Mu is solar units. This is obtained from the halo mass function (dashed lines) by adding the global subhalo mass 
function and subtracting the mass function of multiply-occupied halos (or equivalently multiplying by the probability of single occupation). 
More details are given in App. A. At z = 0 we also report as a dotted line the resulting cluster and group halo mass function (obtained 
by subtraction of the solid from the dashed line); this is compared with the determinations by Boehringer et al. (2014; circles) from X-ray 
observations of groups and clusters and by Martinez et al. (2002; stars) from optical observations of loose groups. 
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Figure 15. Relationship between BH mass Mbh and halo mass Mh from the abundance matching technique, at redshift z = 0 (orange), 
1 (red), 3 (green), and 6 (blue). Top panel shows the results when a one-to-one (i.e., no scatter) relationship Mbh vs. Mh is assumed, 
while bottom panel shows the resulting average relationship when a Gaussian distribution in Mbh at given Mh with a scatter of 0.4 dex 
(see text) is assumed. In both panels the black errorbar illustrates the typical associated uncertainty, and the dotted lines highlight the 
ranges not covered by the current data on the BH mass function. The grey shaded areas show the relations at z = 0 from the BH mass 
functions uncorrected (dark grey) and corrected (light grey) for pseudobulges by Shankar et al. (2009) and (2012), respectively. 
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Figure 16. Relationship between bolometric AGN luminosity Lagn and halo mass Mh from the abundance matching technique, at 
redshift 2 : = 0 (orange), 1 (red), 3 (green), and 6 (blue). Top panel shows the results when a one-to-one (i.e., no scatter) relationship Lagn 
vs. Mh is assumed, while bottom panel shows the resulting average relationship when a Gaussian distribution in Lagn at given Mh with 
a scatter of 0.4 dex (see text) is assumed; in both panels the black errorbar illustrates the typical associated uncertainty, and the dotted 
lines highlight the ranges not covered by the current data on the AGN luminosity function. 
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Figure 17. The AGN bias as a function of redshift z. Results from the abundance matching technique are illustrated by magenta (Tagn- 
average bias) and yellow (MBH -avera g e bias) continuous lines, with the hatched areas showing the associated uncertainty; specifically, the 
magenta curves refer to different AGN luminosities Lagn > 10 10-5 , 10 12 , and 10 13 ' 5 Lq and the yellow curves to different BH masses 
Mbh > 10 6 , 10 8 , and 10 10 Mq as labeled. Black dotted lines illustrate for comparison the halo bias referring to different halo masses from 
10 9 to 10 13 Mq as labeled. The inset shows the AGN bias from the abundance matching technique at redshift 2: = 2 as a function of the 
bolometric AGN luminosity Lagn- Optical data are from Shen et al. (2009; blue circles), Ross et al. (2009; blue diamonds), Da Angela 
et al. (2008; blue reversed triangles), Myers et al. (2007; blue pentagons), Porciani &; Norberg (2006; blue stars), Croom et al. (2005; 
blue squares), White et al. (2012; blue triangles); X-ray data from Allevato et al. (2011; red triangles), Allevato et al. (2014; red circles), 
Mountrichas et al. (2013; red squares), and Starikova et al. (2012; red diamonds). 
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M«or-M„ Relationship 
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Figure 18. Relationship between stellar mass M s tar and halo mass Mh from the abundance matching technique, at redshift z = 0 
(orange), 1 (red), 3 (green), and 6 (blue). The results refer to the average relationship when a Gaussian distribution in M* at given 
Mh with a scatter of 0.15 dex is assumed (the one-to-one relationship is practically identical); the black errorbar illustrates the typical 
associated uncertainty, and the dotted lines highlight the ranges not covered by the current data on the stellar mass function. The grey 
shaded area shows the relation at z = 0 obtained from the observed stellar mass function by Bernardi et al. (2013). The inset illustrates 
the efficiency M*//b Mh for the conversion of the initial baryonic mass /b Mh = 0.2 Mh associated to the halo into the final stellar mass 
M*. 
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Relationship @ z = 0: comparison 
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Figure 19. Relationship between stellar mass M* and halo mass Mh from the abundance matching technique at redshift z = 0 (orange 
line), when the galaxy halo mass function (dashed) or the full halo mass function (solid) are adopted. The grey shaded areas show the 
relations at 2 : = 0 obtained from the observed stellar mass function by Bernardi et al. (2013), matched with the galaxy (light grey) or 
the overall (dark grey) halo mass function. The green solid line refers to the result by Behroozi et al. (2013), the blue solid line to that 
by Moster et al. (2013), and the magenta dashed line to that by Shankar et al. (2006). Data from gravitational lensing measurements in 
groups and clusters of galaxies are from Han et al. (2014; filled stars), Velander et al. (2014; filled circles), and Mandelbaum et al. (2006; 
filled squares). 
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L SF r (SFR)-M h Relationship 
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Figure 20. Average relationship between bolometric SFR-luminosity Lsfr and halo mass Mh from the abundance matching technique, 
at redshift z = 0 (orange), 1 (red), 3 (green), 6 (blue), and 10 (cyan). A Gaussian distribution in Lsfr at given Mh with a scatter of 
0.15 dex (see text) is assumed (the one-to-one relationship is identical); the black errorbar illustrates the typical associated uncertainty, 
and the dotted lines highlight the ranges not covered by the current data on the SFR-luminosity function. Solid lines refer to the overall 
SFR-luminosity function, while dot-dashed lines to dust-corrected UV luminosity function only. At z = 10 the cyan shaded area for small 
halo masses illustrates the systematic uncertainty related to the faint-end slope of the SFR-luminosity function. 
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Galaxy bias vs. redshift 
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Figure 21. The galaxy bias as a function of redshift z. Results from the abundance matching technique are illustrated by magenta 
(LsFR _avera g e bias) and yellow (M*-average bias) continuous lines, with the hatched areas showing the associated uncertainty; specifically, 
the magenta curves refer to different SFRs> 0.03, 3, and 300 Mq yr —1 and the yellow curves to different stellar masses M* > 10 7 , 10 9 , 
and 10 11 Mq as labeled. Black dotted lines illustrate for comparison the halo bias referring to different halo masses from 10 9 to 10 13 Mq 
as labeled. Data for FIR/sub-mm bright galaxies (filled orange stars) are from Webb et al. (2003), Blain et al. (2004), Weiss et al. (2009), 
Hickox et al. (2012), Bianchini et al. (2014), for FIR/sub-mm faint galaxies (orange empty stars) are from Ono et al. (2014), for passive 
BzK galaxies (green filled triangles) are from Grazian et al. (2006), Quadri et al. (2007), Blanc et al. (2008), Furusawa et al. (2011), 
Lin et al. (2012), for starforming BzK galaxies (green empty triangles) are from Hayashi et al. (2007), Blanc et al. (2008), Furusawa et 
al. (2011), for bright Lyman Break Galaxies (blue filled circles) are from Ouchi et al. (2004), Adelberger et al. (2005), Lee et al. (2006), 
Overzier et al. (2006), for faint Lyman Break Galaxies (blue empty circles) are from Bielby et al. (2013), Barone-Nugent et al. (2014), 
for Lyman-a Emitters (cyan diamonds) are from Gawiser et al. (2007), Ouchi et al. (2010), Guaita et al. (2010), for passively-evolving 
Early-Type Galaxies (red filled squares) are from Hawkins et al. (2003), Guzzo et al. (2008), Georgakakis et al. (2014), for Luminous Red 
Galaxies (red empty squares) are from Tegmark et al. (2006), Ross et al. (2007). 
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Cosmic sSFR 
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Figure 22. Cosmic sSFR as a function of redshift. Solid line refers to the overall SFR-luminosity function, while dot-dashed line to 
dust-corrected UV luminosity function only, and the dotted line illustrates the model by Madau & Dickinson (2014). IR data are from 
Ilbert et al. (2014; red circles, referring to M* > 10 10-5 Mq), Damen et al. (2009; orange squares), Reddy et al. (2008; orange stars), 
Noeske et al. (2007; orange pentagons); UV data are from Gonzalez et al. (2014; cyan circles), Salmon et al. (2015; cyan squares), Stark 
et al. (2013; cyan stars), Feulner et al. (2005; cyan pentagons); radio data are from Karim et al. (2011; yellow circles), Daddi et al. (2007; 
yellow squares). 
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SFR-M SUJf Relotionship 
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Figure 23. Relationship between the SFR and the final stellar mass from the abundance matching technique (the so called ‘main 
sequence’ of starforming galaxies), at redshift z = 0 (orange), 1 (red), 3 (green), and 6 (blue). Results with and without scatter are almost 
undistinguishable. Solid lines refer to the overall SFR-luminosity function, while dot-dashed lines to dust-corrected UV luminosity function 
only. The black errorbar illustrates the typical associated uncertainty. The dotted lines highlight the ranges not covered by the current data 
on the SFR-luminosity and stellar mass functions, or where the determination from the abundance matching technique is largely uncertain 
because of the flatness at the faint end of the stellar mass function. Observational estimates are in the range z ~ 0 — 3 by Speagle et al. 
(2014; orange, red and green areas), at 2 ~ 0 by Renzini & Peng (2015; yellow area), at z ~ 2 by Rodighiero et al. (2014; light grey area), 
and at 2 ~ 6 from UV determinations by Salmon et al. (2015; blue open circles). 







50 


R. AVERSA ET AL. 


M BH — M slor Relationship 
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Figure 24. Relationship between the BH mass Mbh and the stellar mass M* from the abundance matching technique, at redshift z = 0 
(orange), 1 (red), 3 (green), and 6 (blue). Top panel shows the results when a one-to-one (i.e., no scatter) relationship Mbh vs. Mh is 
assumed, while bottom panel shows the resulting average relationship when a Gaussian distribution in Mbh at given Mh with a scatter of 
0.3 dex (see text) is assumed; in both panels the black errorbar illustrates the typical associated uncertainty, and the dotted lines highlight 
the ranges not covered by the current data on the BH and stellar mass functions. The shaded areas show the relations at z = 0 obtained 
from the matching between the stellar and the BH mass functions, uncorrected (dark grey) and corrected (light grey) for pseudobulges. 
Data points are from the compilation by Haring & Rix (2004), with the dashed line representing their best-fit relation; the relation proposed 
by Kormendy & Ho (2013) is also shown as a dot-dashed line. The inset illustrates the corresponding BH-to-stellar mass ratio Mbh/M* 
as a function of M*. 
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Cosmic mass densities 
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Figure 25. The evolution with redshift of the mass density in the overall DM halo population (dotted green line), in galactic DM halos 
(solid green line), in stars (solid orange line), and in black holes (solid cyan line). In the inset the luminosity density ratio Pl sfr /pl agn 
(blue line) and the mass density ratio pM+ /Pm bh (red line) are illustrated. 
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Subholo Mass Function 
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Figure Al. The subhalo mass function JV(logm) vs. the ratio between the satellite and the halo masses m/Mn, computed according to 
the prescriptions by van den Bosch & Jiang (2014). The black line refers to the unevolved mass function, and colored lines to the evolved 
mass function at 2 : = 0 (orange), 1 (red), 3 (green), and 6 (blue). At 2 : = 0 the solid line refers to a mass of the host of Mh = 10 13 Mq, 
dashed line to a 10 12 Mq, and dotted to 10 14 Mq. 
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Holo & Subhalo Mass Function 
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Figure A2. The overall contribution of subhalos to the halo mass function, vs. the halo mass Mh- Solid lines show the halo mass function 
including subhalos, dashed lines show the halo mass function without subhalos, and dotted lines show the subhalo contribution. Results 
are plotted at redshift z = 0 (orange), 1 (red), 3 (green), and 6 (blue). 
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Halo Occupation Number 
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Figure A3. The halo occupation number vs. the host halo mass Mh- Results are plotted at redshift z = 0 (orange), 1 (red), 3 (green 
and 6 (blue). At 2 : = 0 solid line refers to a minimum satellite mass m m i n = 10 11 Mq, dashed to 10 10 ' 5 Mq, and dotted to 10 11 ' 5 Mq. 
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SMBH Mass Function: effect of dry merging 
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Figure Bl. Effect of dry mergers on the late z < 1 evolution of the supermassive BH mass function. The red line represents the mass 
function at redshift z ~ 1 (at higher redshift dry merging effects are negligible), while the other colored lines illustrates its evolution toward 
z ~ 0 due to merging and in-situ accretion. Specifically, the BH merging rate is assumed to mirror the DM merging rates as given by 
Stewart et al. (2009) for major mergers (yellow line) and for minor mergers (pink line); the result for in-situ accretion only is plotted for 
reference as an orange line. Data points and shaded areas as in Fig. [4] 
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Stellar Mass Function: effect of merging 
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Figure B2. Effect of dry mergers on the late z < 1 evolution of the galaxy stellar mass function as derived from the continuity equation. 
The red line represents the mass function at redshift z ~ 1 (at higher redshift dry merging effects are negligible), while the other colored 
lines illustrate its evolution toward z ~ 0 due to merging and in-situ formation. Specifically, the galaxy merging rate is computed according 
to Stewart et al. (2009) for major mergers (yellow line) and for minor mer gers (pink line); the result for in-situ formation only is plotted 
for reference as an orange line. Data points and shaded areas as in Fig. 1101 
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UV Luminosity Function 
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Figure Cl. The SFR-luminosity function 7V(log Lsfr) at redshift z = 0 (orange), 1 (red), 3 (green), and 6 (blue), vs. the luminosity 
Lsfr associated to the SFR (lower axis) and vs. the SFR (upper axis). Solid lines is our rendition of the luminosity function based on 
the UV data at the faint and FIR data at the bright; this is the same plotted in Fig. |9] Dot-dashed lines is a rendition based only on 
dust-corrected UV data, and dashed lines from dust-uncorrected UV data. The inset shows the corresponding SFR-luminosity densities. 
Data points have been omitted for clarity. 
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Stellar Mass Function 



log M stor [M e ] 


Figure C2. Effect of adopting the SFR-luminosity functions obtained basing on FIR and UV data (solid lines), only dust-corrected UV 
data (dot-dashed lines) and only dust-uncorrected UV data (da shed lines) as input of the continuity equation to obtain the stellar mass 
function. The solid lines are the same outputs shown in Fig. 1101 to be in very good agreement with data points (here omitted for clarity). 








Table 1 

Input/Output Generalized Schechter Functions. 


function 

log < 1-0 

k& 1 

k$ 2 

fc# 3 

log Vq 

kxi 

kx 1 

kx 3 

a(z 0 ) 

ka l 

k a 2 

k a 3 

w(«o) 

kui 

ku, 2 

ku, 2 


AGN LF 

-3.80 

0.45 

-1.00 

0.00 

10.90 

1.10 

6.94 

-11.55 

1.40 

-1.70 

3.40 

-1.75 

0.36 

0.62 

-1.59 

0.8 

SFR LF 

-2.40 

-2.30 

6.20 

-4.90 

10.90 

3.20 

-1.40 

-2.10 

1.20 

0.50 

-0.50 

0.20 

0.70 

-0.15 

0.16 

0.00 

BH MF 

-2.30 

-0.40 

-1.80 

-1.50 

8.07 

-0.80 

7.30 

-9.20 

1.35 

-0.10 

0.40 

0.30 

0.46 

0.05 

0.18 

-0.55 

Stellar MF 

-2.10 

-0.80 

1.65 

-3.10 

10.85 

0.00 

0.00 

-1.90 

1.20 

0.00 

-0.40 

0.55 

0.65 

0.00 

-0.40 

0.55 

HMF 

-3.97 

0.00 

0.00 

1.50 

14.00 

-0.90 

-1.90 

-1.10 

1.80 

0.50 

0.10 

0.70 

0.47 

0.45 

-0.10 

-0.45 

GHMF 

-3.35 

0.50 

0.1 

-1.50 

13.05 

-0.80 

0.00 

-1.30 

1.88 

0.30 

-0.40 

1.30 

1.10 

-0.10 

0.00 

-0.43 


Note. — Typical tolerance on the parameters is less than 10%. See Sect. |2.1.ll for details on the redshift evolution of the parameters. 
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Table 2 

Fits to Abundance Matching Results 


function 

log ATq 

fcjvi 

fcjV2 

fejV3 

log Mb o 

kMi 

kM2 

*M3 

<20 

kczi 

k a 2 

foa 3 

UJQ 


ku,2 

kcj3 







Y = 

N x [ 

;*r+( 

1 

3 









Lagn - Mu (ltol) 

12.25 

-1.90 

1.45 

1.10 

12.30 

0.00 

0.00 

0.00 

-1.65 

1.20 

-4.50 

0.60 

-0.65 

-3.50 

-0.40 

0.80 

Lagn ~ M h (aver.) 

12.33 

-2.00 

2.70 

-1.30 

12.50 

0.00 

0.00 

0.00 

-1.60 

0.00 

-3.30 

1.30 

-0.50 

-2.10 

-1.60 

2.60 

Lsfb. — Mu (aver.) 

11.25 

0.00 

0.00 

0.00 

12.20 

-1.20 

0.00 

0.00 

-1.30 

-3.00 

-0.50 

1.20 

-0.50 

-1.50 

0.00 

1.50 

A/bh — Mu (ltol) 

8.00 

-0.40 

0.70 

-0.80 

11.90 

0.00 

0.00 

0.00 

-1.10 

-0.80 

-1.50 

0.10 

-1.10 

-0.80 

-1.50 

0.10 

Mbh — Mr (aver.) 

8.20 

-0.20 

0.80 

-1.50 

12.20 

0.00 

0.00 

0.00 

-1.40 

-1.30 

-0.30 

0.10 

-0.80 

-0.40 

-1.10 

0.10 

M* — Mr (aver.) 

10.40 

-0.80 

0.80 

-0.20 

11.50 

0.00 

0.00 

0.00 

-2.20 

-1.90 

-1.60 

4.70 

-0.75 

-0.30 

-1.80 

2.60 

SFR-M* (aver.) 

1.90 

0.00 

0.00 

0.00 

11.60 

-1.90 

-2.50 

1.70 

-1.60 

0.00 

1.50 

-0.20 

-0.50 

-1.70 

3.50 

-2.00 






Y 

= N x 

[(&)■ + 

3 









Mbh — M* (ltol) 

7.33 

0.10 

-0.70 

1.00 

10.60 

0.00 

0.00 

0.00 

1.60 

0.00 

-0.20 

2.80 

0.60 

0.40 

-0.40 

2.20 

M B u - M* (aver.) 

7.15 

0.00 

-0.60 

0.00 

10.50 

0.00 

0.00 

0.00 

1.30 

0.30 

0.00 

0.90 

0.60 

0.20 

0.40 

0.90 


Note. — Typical tolerance on the parameters is less than 10%. See Appendix D for details on the redshift evolution of the parameters. 
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Table 3 

Main Results 


Results Sections Figures 


AGN luminosity function 
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BH mass function 
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SFR luminosity function 

22T] 

19] 

Stellar mass function 
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Galaxy halo mass function 
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AGN/BH bias 
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M* — Mh relationship 
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Galaxy bias 
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